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ABSTRACT 

We present a comprehensive set of convergence tests which explore the role of various 
numerical parameters on the equilibrium structure of a simulated dark matter halo. 
■ We report results obtained with two independent, state-of-the-art, multi-stepping, 

t^J- \ parallel N-body codes: PKDGRAV and GADGET. We find that convergent mass profiles 

l/^ . can be obtained for suitable choices of the gravitational softening, timestep, force ac- 

curacy, initial redshift, and particle number. For softenings chosen so that particle 
C> discreteness effects are negligible, convergence in the circular velocity is obtained at 

04 . radii where the following conditions are satisfied: (i) the timestep is much shorter than 

the local orbital timescale; (ii) accelerations do not exceed a characteristic acceleration 
. imprinted by the gravitational softening; and (iii) enough particles are enclosed so that 

O ^ the collisional relaxation timescale is longer than the age of the universe. Convergence 

also requires sufficiently high initial redshift and accurate force computations. Poor 
spatial, time, or force resolution leads generally to systems with artificially low central 
density, but may also result in the formation of artificially dense central cusps. We have 
\ explored several adaptive time-stepping choices and obtained best results when indi- 

vidual timesteps are chosen according to the local acceleration and the gravitational 
softening (Ati oc (e/ai) 1 / 2 ), although further experimentation may yield better and 
more efficient criteria. The most stringent requirement for convergence is typically that 
; ^ ■ imposed on the particle number by the collisional relaxation criterion, which implies 

that in order to estimate accurate circular velocities at radii where the density contrast 
may reach ~ 10 6 , the region must enclose of order 3000 particles (or more than a few 
times 10 6 within the virial radius). Applying these criteria to a galaxy-sized ACDM 
halo, we find that the spherically-averaged density profile becomes progressively shal- 
lower from the virial radius inwards, reaching a logarithmic slope shallower than —1.2 
at the innermost resolved point, r ~ 0.005 t^oo, with little evidence for convergence to 
a power-law behaviour in the inner regions. 
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1 INTRODUCTION 

Over the past few decades, cosmological N-body simulations 
have led to impressive strides in our understanding of struc- 
ture formation in universes dominated by collisionless dark 
matter. Such simulations have provided an ideal test-bed 
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for analytic theories of structure formation, and have been 
used to validate and motivate a variety of theoretical in- 
sights into the statistics of hierarchical clustering (e.g., Press 
& Schechter 1974, Bardeen et al. 1986, Bond et al. 1991, 
Lacey & Cole 1993, Mo & White 1996). In particular, N- 
body simulations have played a pivotal role in providing a 
clear framework within which the CDM cosmogony may be 
compared with observation, and in establishing Cold Dark 
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Matter (CDM) as the leading theory of structure formation 
(Davis et al. 1985). 

This work has led to the development of a robust theo- 
retical framework which provides an accurate statistical de- 
scription of structure growth through gravitational instabil- 
ity seeded by Gaussian primordial density fluctuations. It 
is now possible to predict with great accuracy, and based 
only on the initial power spectrum of the primordial fluc- 
tuations, a number of important statistics that characterize 
the large scale structure of the universe; e.g., the mass func- 
tion and clustering of dark matter halos and their evolution 
with redshift (e.g., Jing 1998, Sheth & Tormen 1999, Jenk- 
ins et al. 2001) the non-linear evolution of the dark matter 
power spectrum and correlation functions (e.g., Hamilton et 
al. 1991, Peacock & Dodds 1996), as well as the topological 
properties of the large scale structure (e.g., Gott, Weinberg 
& Melott 1987). 

The impact of such simulation work has been greatest 
in the non-linear regime, where analytic calculations offer 
little guidance. Recently, and as a result of the development 
of efficient algorithms and of the advent of massively paral- 
lel computers, it has been possible to apply N-body studies 
to the investigation of structure on small, highly non-linear 
scales. These studies can now probe scales comparable to 
the luminous radii of individual galaxies, thus enabling di- 
rect comparison between theory and observation in regions 
where luminous dynamical tracers are abundant and easi- 
est to observe. Predicting the structure of dark matter halos 
on kpc and sub-kpc scales, where it can be compared di- 
rectly with observations of galactic dynamics, is one of the 
premier goals of N-body experiments, and there has been 
steady progress in this area over the past few years. 

Building upon the early work of Frenk et al. (1985, 
1988), Quinn, Salmon & Zurek (1986), Dubinski & Carl- 
berg (1991) and Crone, Evrard & Richstone (1993), Navarro, 
Frenk & White (1996, 1997, hereafter NFW) found that, 
independently of mass and of the value of the cosmolog- 
ical parameters, the density profiles of dark matter halos 
formed in various hierarchical clustering cosmogonies were 
strikingly similar. This 'universal' structure can be charac- 
terized by a spherically-averaged density profile which differs 
substantially from the simple power law, p(r) oc r~P, pre- 
dicted by early theoretical studies (Gunn & Gott 1972, Fill- 
more & Goldreich 1984, Hoffmann & Shaham 1985, White & 
Zaritsky 1992). The profile steepens monotonically with ra- 
dius, with logarithmic slopes shallower than isothermal (i.e. 
(3 < 2) near the centre, but steeper than isothermal (/3 > 2) 
in the outer regions. 

NFW proposed a simple formula, 

p( r ) = 5c m 

Pent (r/r s )(l + r/r a y { ' 

which describes the density profile of any halo with only 
two parameters, a characteristic density contrast^, S c , and 
a scale radius, r s . Defining the mass of a halo as that con- 
tained within r2oo, the radius of a sphere of mean density 

t We use the term 'density contrast' to denote densities expressed 
in units of the critical density for closure, p cr it = 3H 2 /8nG. We 
express the present value of Hubble's constant as H(z = 0) = 
H = 100 h km s" 1 Mpc" 1 



contrast 200, there is a single adjustable parameter that fully 
describes the mass profile of halos of given mass: the 'con- 
centration' ratio c = r2oo/r 3 - 

For the sake of this discussion, the two main points to 
note from the work of NFW are the following: (i) the density 
profile in the inner regions of the halo is shallower, and in the 
outer regions steeper, than isothermal, and (ii) there is no 
well defined value for the central density of the dark matter, 
which can in principle climb to arbitrarily large values near 
the centre. 

Conclusion (i) is important, since it is a feature of dark 
halo models that is required by observations. For example, 
it implies that the characteristic speeds of dynamical trac- 
ers may be lower near the centre than in the main body of 
the system, as observed in disk galaxies, where the veloc- 
ity dispersion of the bulge is lower than indicated by the 
maximum rotation speed of the surrounding disk, as well 
as in galaxy clusters, where the velocity dispersion of stars 
in the central cluster galaxy is lower than that of the clus- 
ter as a whole. Conclusion (ii) is also important, since there 
have been a number of reports in the literature arguing that 
the shape of the rotation curves of many disk galaxies rules 
out steeply divergent dark matter density profiles (Flores & 
Primack 1994, Moore 1994, de Blok et al. 2001, but see van 
den Bosch & Swaters 2001), a result that may signal a gen- 
uine crisis for the CDM paradigm on small scales (see, e.g., 
Sellwood & Kosowsky 2000, Moore 2001). 

These general results of the work by NFW have been 
confirmed by a number of subsequent studies (Cole & Lacey 
1996, Fukushige & Makino 1997, Huss, Jain & Steinmetz 
1999, Moore et al. 1998, Jing & Suto 2000), although there 
is some disagreement about the innermost value of the loga- 
rithmic slope. Moore et al. (1998), Ghigna et al.(2000), and 
Fukushige & Makino (1997, 2001) have argued that den- 
sity profiles diverge near the centre with logarithmic slopes 
considerably steeper than the asymptotic value of (3 = 1 in 
NFW's formula. Kravtsov et al. (1998), on the other hand, 
initially obtained much shallower inner slopes (/3 ~ 0.7) 
in their numerical simulations, but have now revised their 
conclusions; these authors now argue that CDM halos have 
steeply divergent density profiles but, depending on evolu- 
tionary details, the slope of a galaxy-sized halo at the in- 
nermost resolved radius may vary between —1.0 and —1.5 
(Klypin et al. 2001). 

Since steep inner slopes are apparently disfavoured by 
rotation curve data it is important to establish this re- 
sult conclusively; if confirmed, it may offer a way to fal- 
sify the CDM paradigm on small scales. Unfortunately, ob- 
servational constraints are strongest just where theoretical 
predictions are least trustworthy. For example, the alleged 
disagreement between observed rotation curves and cuspy 
dark halo models is most evident for sub-L* galaxies on 
scales of ~ 1 h^ 1 kpc or less. For typical circular speeds of 
~ 100 km s" 1 , this corresponds to regions where the den- 
sity contrast exceeds ~ 10 . Orbital times in these regions 
are of order 10 of the age of the universe, implying that 
N-body codes must be able to follow particles accurately for 
several thousand orbits. Few cosmological codes have been 
tested in a systematic way under such circumstances. Fur- 
thermore, the cold dark matter halos that host typical disk 
galaxies are thought to extend out to a few hundred kpc, im- 
plying that the ~kpc scale probed by observations involves a 
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very small fraction of the mass and volume of the dark halo. 
As a consequence, these regions are vulnerable to numeri- 
cal artifacts in N-body simulations stemming, for example, 
from the gravitational softening or the number of particles. 

Extreme care is thus needed to separate numerical ar- 
tifacts from the true predictions of the Cold Dark Matter 
model. In order to validate or 'rule out' the CDM cosmogony 
one must be certain that model predictions on the relevant 
scales are accurate, robust, and free of systematic numer- 
ical uncertainties. Although there have been some recent 
attempts at unravelling the role of numerical parameters on 
the structure of simulated dark matter halos, notably in the 
work of Moore et al. (1998), Knebe et al. (2000), Klypin et 
al. (2001) and Ghigna et al. (2000), the conclusions from 
these works are still preliminary and, in some cases, even 
contradictory. 

To cite an example, Moore et al. (1998) argue that the 
smallest resolved scales correspond to about half the mean 
inter-particle separation within the virial radius, and con- 
clude that many thousands of particles are needed to re- 
solve the inner density profile of dark matter halos. Klypin 
et al. (2001), on the other hand, conclude that mass profiles 
can always be trusted down to the scale of the innermost 
~ 200 particles, provided that other numerical parameters 
are chosen wisely. Ghigna et al. (2000) suggest an additional 
convergence criterion based on the gravitational softening 
length scale, and argue that convergence is only achieved 
on scales that contain many particles and that are larger 
than about ~ 3 times the scale where pairwise forces be- 
come Newtonian. Understanding the origin of such disparate 
conclusions and the precise role of numerical parameters is 
clearly needed before a firm theoretical prediction for the 
structure of CDM halos on ~kpc scales may emerge. 

Motivated by this, we have undertaken a large series of 
numerical simulations designed to clarify the role of numeri- 
cal parameters on the structure of simulated cold dark mat- 
ter halos. In particular, we would like to answer the following 
question: what regions of a simulated dark matter halo in 
virial equilibrium can be considered reliably resolved? This 
question is particularly difficult because of the lack of a the- 
ory with which the true structure of dark halos may be pre- 
dicted analytically, so the best we can do is to establish the 
conditions under which the structure of a simulated dark 
halo is independent of numerical parameters. This is the 
question which we endeavor to answer in this paper. 

There is a long list of considerations and numerical pa- 
rameters that may influence the structure of simulated dark 
halos: 

• the N-body code itself 

• the procedure for generating initial conditions 

• the accuracy of the force computation 

• the integration scheme 

• the initial redshift 

• the time-stepping choice 

• the gravitational softening 

• the particle number 

Clearly the list could be substantially longer, but the items 
above are widely considered the most important concern- 
ing the structure of simulated dark halos. Before we pro- 
ceed to analyze their role, we must decide which properties 
of a dark matter halo we will assess for numerical conver- 



gence. Because, as mentioned above, disk galaxy rotation 
curves seem to pose at present one of the most pressing 
challenges to the CDM paradigm on small scales, we have 
decided to concentrate on the spherically-averaged mass pro- 
file, as measured b y the radi al dependence of the circular 
velocity, V c (r) = GM(r)/r, or, equivalently, by the inner 
mean density profile, p(r) = 3M(r)/4-7rr 3 . 

We note that the convergence criteria derived here ap- 
ply strictly only to these properties, and that others, such as 
the three-dimensional shape of halos, their detailed orbital 
structure, or the mass function of substructure halos, may 
require different convergence criteria. 

The basic philosophy of our convergence testing proce- 
dure is to select a small sample of halos from a cosmological 
simulation of a large periodic box and to resimulate them 
varying systematically the parameters listed above, search- 
ing for regions in parameter space where the circular velocity 
curves are independent of the value of the numerical param- 
eters, down to the smallest scales where Poisson uncertain- 
ties become important, i.e., roughly down to the radius that 
contains ~ 100 particles. 

Overall, this is a fairly technical paper of interest mostly 
to practitioners of cosmological N-body simulations. Readers 
less interested in numerical details may wish to skip to § 6, 
where we discuss in detail the converged inner mass profile of 
the galaxy-sized ACDM halo used in our convergence study. 
The more technical sections include: 

• a discussion of the N-body codes used in this work, ini- 
tial conditions setup, and analysis procedure (§ 2 and Ap- 
pendix) 

• a general discussion of the consequences of discrete- 
ness effects on simulations of dark matter halos, including a 
derivation of "optimal" choices (for given particle number) 
of the timestep and the gravitational softening (§3) 

• a comparison between single- and multi-timestepping 
techniques (§ 4) 

• a discussion of the role of the gravitational softening, 
the initial redshift, the force accuracy, and the particle num- 
ber on the inner mass profile of simulated halos (§ 5) 

Finally, a worked example of how to choose optimal param- 
eters for a high-resolution simulation is presented in § 5.5. 
We summarize our main conclusions in § 7. 



2 NUMERICAL METHODS 
2.1 N-body Codes 

Most simulations reported in this paper have been per- 
formed with the parallel N-body code GADGET, writ- 
ten by Volker Springel, and available from http:// 
www.mpa-garching.mpg.de/gadget (Springel, Yoshida & 
White 2001). In order to test the dependence of our results 
on the particular algorithmic choices made in GADGET, we 
have also used PKDGRAV, a code written by Joachim Stadel 
and Thomas Quinn (Stadel 2001). As we discuss in § 3 and 
§ 5, the two codes give approximately the same results for 
appropriate choices of the numerical parameters. We have 
not attempted to carry out a detailed comparison of the 
relative efficiency or speed of the codes; such comparison 
is heavily dependent on the particular architecture of the 
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hardware used, and on a variety of optimization and tuning 
procedures. We do note, however, that neither code seems 
obviously to outperform the other when our strict numerical 
convergence criteria are met. 

The two N-body codes share a number of similarities. 
They both evaluate accelerations ('forces') on individual 
particles due to all others using a hierarchical tree data 
structure (Barnes & Hut 1986, Jernigan & Porter 1989), 
and (optionally) use individually adaptive time-stepping 
schemes to advance the integration of each particle. Periodic 
boundary conditions are handled in both codes via Ewald's 
summation technique (Hernquist, Bouchet & Suto 1991), al- 
though the implementation of the algorithm in each code is 
different . 

Gravitational softening is introduced in the form of a 
'spline' mass distribution (see, e.g., Hernquist & Katz 1989, 
Navarro & White 1993) which, unlike the more traditional 
'Plummer' softening of the early generation of N-body codes 
(see, e.g., Aarseth 1985), converges for pairwise interactions 
exactly to the Newtonian regime at a finite radius. The 
length scale of the spline kernel, e», can be chosen individu- 
ally for each particle in PKDGRAV. GADGET, on the other hand, 
allows for different softenings to be chosen for up to six dif- 
ferent particle 'species'. We quote the values of e; so that 
gravitational interactions between two particles are fully 
Newtonian for separations larger than 2e;. 

The codes differ substantially in their implementation 
of the tree construction, in the force-evaluation algorithms 
and in the integrator scheme. Whereas PKDGRAV uses a spa- 
tial binary tree for gravity calculations, GADGET uses a ver- 
sion of the Barnes-Hut geometric oct-tree. Distant tree- 
node contributions to the force calculations include up to 
quadrupole expansion terms in GADGET, but up to hexade- 
capole in PKDGRAV. The tree is rebuilt every timestep in the 
version of PKDGRAV that we tested (although this is not the 
case in the most up-to-date version), whereas we rebuild 
the tree in GADGET dynamically after ~ 0.1 Ntot force com- 
putations since the last full reconstruction. (N to t is the total 
number of particles in the simulation.) 

Finally, GADGET uses a simple second-order DKD (drift- 
kick-drift) leap-frog integrator scheme with expansion fac- 
tor as the integration variable, whereas PKDGRAV adopts a 
cosmic time-based KDK (kick-drift-kick) algorithm. All in- 
tegrations are carried out in comoving coordinates. Details 
of these codes may be found in Springel et al. (2001), and 
in Stadel (2001). In the following subsections we describe 
the numerical setup used for the two codes. All simulations 
have been run on the IBM-SP supercomputer facilities at 
the University of Victoria (Canada), and on the T3Es at 
the Edinburgh Parallel Computer Centre (U.K.) and at the 
Max-Planck Rechenzentrum in Garching (Germany). 

2.1.1 GADGET 

GADGET has been the main simulation code used in this study, 
and it evolved as the project unfolded from the first public 
release vl.O to the latest available release vl.l. All of the 
results presented here have been obtained with the latest 
version of the code. 

GADGET presents the user with a number of options re- 
garding time-stepping choices and the accuracy of the force 
calculations. In all cases we have used the tree node-opening 



criterion recommended by Springel et al. (2001), where a 
Barnes-Hut opening criterion with 6 — 0.6 is used for the 
first force computation and a dynamical updating criterion 
is used subsequently. In this criterion, a node is opened if 
MZ 4 >/ acc o-oid t , where M is the mass of the node, I is the 
node-side length, and a id is the acceleration that the par- 
ticle experienced in the previous timestep. The parameter 
/ acc (called ErrTolForceAcc in GADGET'S parameter list) is 
set to 10~ 3 in our standard calculations. This condition can 
be overridden if the -DBMAX compile-time flag is activated. 
Enabling this flag imposes an additional condition for node- 
opening: multipole expansion of a node is only used if, in 
addition to the previous condition, the particle is guaran- 
teed to lie outside the geometric boundaries of the node in 
question^. The results reported in § 5.3 indicate that these 
choices are important to ensure convergence: resolving the 
inner structure of dark halos requires highly accurate forces. 

GADGET uses an integrator with completely flexible 
timesteps. The code carries, for each particle, a time, U, 
position, r», velocity, v,, acceleration, a^, gravitational soft- 
ening, £», and, optionally, a local density, p,, and a local one- 
dimensional velocity dispersion, <7j. From these quantities, 
timesteps, Ati, can be computed for each particle according 
to several possible choices: 



At 



if DtCrit=0 
if DtCrit=l 
if DtCrit=2 
if DtCrit=3 
if DtCrit=4 



(2) 



r)a<y{cTi/ai), 

v P (G Pl y 1/2 , 

f) rTp ram[(Gp i ) 1/2 , (ffi/oj)], 

where DtCrit refers to the runtime input parameter 
ErrTolIntAccuracy in GADGET, and 77 is a dimensionless con- 
stant that controls the size of the timesteps (except for r) a , 
which has dimensions of velocity) §. For ease of reference, we 
shall refer to the various choices for DtCrit using the follow- 
ing mnemonic shorthand: EpsAcc for DtCrit=0; VelAcc for 
DtCrit=l; SgAcc for DtCrit=2; SqrtRho for DtCrit=3; and 
RhoSgAcc for DtCrit=4, respectively 

We report below results obtained with several of these 
choices. Unless specified, a maximum timestep was imposed 
so that all particles took at least 200 timesteps during the 
whole integration. In practice, this limit affects a very small 
fraction of the particles in a typical run: resolving the inner 
structure of dark halos requires typically several thousand 
timesteps. 

2.1.2 PKDGRAV 

In the PKDGRAV runs reported below we have only explored 
variations in two parameters: the time-stepping parameter, 
77, and the gravitational softening, ti. We note, however, that 
PKDGRAV is a very flexible code that includes a number of 
choices for the integrator scheme and time-stepping, and we 
have by no means explored all of its options. PKDGRAV was 
mainly used in this study to verify that the results obtained 
with GADGET are independent of the code utilized. 

t A similar condition is activated by default in PKDGRAV 
§ For convenience we have defined r} a t to be directly proportional 
to the size of the timestep in all cases. For DtCrit=0, n^ e = 
2 X ErrTolIntAccuracy 
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All PKDGRAV simulations that used individual timesteps 
were evolved to z = using 50 system timesteps. The sys- 
tem timestep, AT, is the maximum allowed for any particle. 
Individual particle timesteps are binned in a hierarchy so 
that AU = AT/2", where n was allowed to take any value 
in the range (0,20). This allows particles to take up to ~ 10 s 
timesteps in a run, which means that in practice no signifi- 
cant restrictions have been placed on the minimum timestep. 

Individual particle timesteps were chosen in PKDGRAV 
runs in a manner analogous to GADGET'S EpsAcc criterion, 
i.e., Ati < rf^jejoi, although quantitatively accelerations 
differ because of the choice of integration variables. The 
parameter rj specifies the size of the timesteps and, conse- 
quently, the overall time accuracy of the integration. Finally, 
the force accuracy in PKDGRAV is controlled by 9, a redshift- 
dependent opening-node criterion. We have chosen for all 
runs 9 = 0.55 (2 > 2) and 9 = 0.7 for 2 < 2. 



2.2 The Initial Conditions 

Setting up initial conditions that faithfully represent the cos- 
mogony one wishes to investigate is a crucial step in the 
simulation process and, despite the popularity of cosmolog- 
ical N-body simulations, there is surprisingly little detail 
in the literature regarding how this is tackled by different 
groups. The major references on this topic in the refereed 
literature are the work of Efstathiou et al. (1985) and the re- 
cent papers by Bertschinger & Gelb (1991), Pen (1997), and 
by Bertschinger (2001; see also http://arcturus.mit.edu/ 
cosmics and http://arcturus.mit.edu/grafic). 

Our particular procedure follows closely that described 
in Efstathiou et al. (1985) and is described in detail in the 
Appendix. It aims to provide a particle realization of a Gaus- 
sian density field with the chosen primordial power spec- 
trum, P(k), on scales and at redshifts where linear theory is 
applicable. 

We adopt the ACDM cosmological model, a low-density 
universe of flat geometry whose dynamics is dominated at 
present by a cosmological constant, flo = 0.3, 57a = 0.7 and 
h = 0.65. We shall assume that the initial power spectrum is 
Harrison-Zel'dovich (P(k) oc k), modified by an appropriate 
cosmological transfer function, T(k). For ACDM simulations 
we have chosen to use the analytic representation of the 
transfer function proposed by Bardeen et al. (1986) with 
shape parameter T = 0.2. 

Our simulations proceed in two stages. Firstly, a large, 
low-resolution, periodic box is run to z — and used to 
select halos targeted for resimulation at much higher resolu- 
tion (consult the Appendix for details). For the first step, we 
have generated a Fourier representation of the fluctuation 
distribution on a 128 3 mesh and have computed displace- 
ments for 128 3 particles initially arranged on a cubic grid. 
The displacements assume an initial redshift of zt = 49 in 
the ACDM cosmogony and are normalized so that at 2 = 
the linear rms amplitude of mass fluctuations on spheres 
of radius 8/i _1 Mpc is as = 0.9. The size of the box is 
ibox = 32.5 ft -1 Mpc (comoving), and the particle mass 
is m„ = 4.55 x lO^o^M©. The dashed curve in Fig- 
ure 1 shows that the power spectrum computed from the 
displaced positions of the 128 3 particles within this box is 




log 10 (k/hM P c-«) 

Figure 1. The dotted line shows the theoretical ACDM power 
spectrum at redshift 2 = 49. The short dashed curve shows the 
measured power spectrum from the initial conditions of the par- 
ent simulation (L box = 32.5 h^ 1 Mpc, JV box = 128 3 ). The solid 
line shows the power spectrum within the high-resolution box se- 
lected for resimulation (L s box = 5.08 ft -1 Mpc, N s ^ ox = 256 3 ). 
The agreement with the theoretical power spectrum is good 
over nearly three orders of magnitude in wavenumber and seven 
decades in amplitude. Significant departures are expected for both 
curves at low k as the number of long-wavelength modes is small. 
The charge assignment scheme causes a small drop at high-fc for 
both curves. The vertical long dashed line marks the scale in the 
rcsimulatcd initial conditions which corresponds to the transition 
between the long waves which are present with the same phase 
and amplitude as the parent simulation and the additional short 
waves added to improve the resolution. See Appendix for more 
details of the computation of the power spectrum. 



in very good agreement with the theoretical power spectrum 
(dotted lines). 

The second stage of the initial conditions generating 
procedure involves selecting a small region within the large 
periodic box destined to collapse into a halo selected for res- 
imulation at higher resolution. In the case we consider here, 
this region is a box of Z/ s box = 5.08 h^ 1 Mpc on a side. The 
advantage of this procedure is that one can in principle in- 
clude many more particles in the high-resolution box than 
were present in the parent simulation (we use AT s b ox = 256 3 
in the case we consider here, giving a highest-resolution par- 
ticle mass of 6.5x 10 5 hT 1 M©). A new Fourier representation 
of the theoretical power spectrum is then generated, retain- 
ing the phases and amplitudes of the Fourier components in 
the parent simulation and adding waves of higher frequency, 
periodic in the high-resolution box, up to the Nyquist fre- 
quency of the high-resolution particle grid. The solid line 
in Figure 1 shows that the power spectrum measured di- 
rectly from particle displacements in the high-resolution box 
is again in good agreement with the theoretical expectation. 

Figure 1 thus demonstrates that the power spectrum is 
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Figure 2. Particle distribution in the initial conditions of our 
128 3 runs at 2; = 49. Clearly seen is the 'amoeba'-shaped region 
containing the highest resolution particles, which is embedded 
within the 5.08 h~ x Mpc high-resolution cube used for the res- 
imulation. Beyond the boundaries of the high-resolution cube lie 
massive particles that coarsely sample the entire volume of the 
32.5 ft -1 Mpc periodic box. 

reproduced well by both the parent simulation and the res- 
imulated region. Altogether, the power spectrum is fit well 
over nearly three decades in wavenumber and seven decades 
in power. The maximum difference between the theoretical 
power spectrum and the measured power spectra is less than 
0.05 dex, except at low wavenumbers where the small num- 
ber of modes makes the variance of the measurement large. 

Outside the high-resolution box, we resample the parti- 
cle distribution in the parent simulation in order to provide 
for the tidal forces which act on the high resolution particles. 
The resampling procedure bins particles into cells whose size 
varies approximately in proportion to their distance from 
the high resolution patch, greatly reducing the total num- 
ber of tidal particles needed to represent the tidal field. Not 
all particles in the high-resolution box will end up near the 
system of interest, so the location on the original grid of se- 
lected particles is used to identify an 'amoeba-shaped' region 
within the cube that is retained at full resolution. Regions 
exterior to the 'amoeba' are coarse sampled with particles 
of mass increasing with distance from the region of interest 
(Figure 2). 

2.3 The Simulations 

The initial conditions file containing the displacement field 
for iV s box = 256 3 particles generated in the way described 
above can be easily rescaled to generate realizations of each 
system with varying particle number or starting redshift. To 
modify the starting redshift, we simply rescale the displace- 
ments and velocities according to the linear growth factor. 



To reduce the particle number, we average successively dis- 
placements in the high-resolution box over 8 neighboring 
cubic cells. We refer to these 'reduced' initial conditions us- 
ing the total number of particles in the high-resolution box: 
256 3 , 128 3 , 64 3 , and 32 3 , respectively (Table 1). 

These realizations may be used to test how numerical 
parameters affect the equilibrium structure of the dark halo 
at z = 0. Since runs with 32 3 particles are relatively inex- 
pensive, we have used them for a large series of simulations 
varying systematically all the numerical parameters under 
scrutiny. This series (which contains several hundred runs) 
allows us to survey the large available parameter space and 
to draw preliminary convergence criteria that are then con- 
firmed with a series of runs with 64 3 particles. The 128 3 
and 256 3 simulations are too expensive to allow a full con- 
vergence study, so fewer of them were carried out, typically 
using values of the numerical parameters close to conver- 
gence. These are used mainly to test the dependence of our 
results on the total number of particles in the simulations. 



2.4 The Halo 

We concentrate our analysis on a single halo selected from 
our sample, although similar runs on two other halos confirm 
the conclusions presented here. The mass accretion history 
of this system is presented in Figure 3. The halo accretes half 
of its present-day mass by z m 0.66 (expansion factor a — 
0.6), when it undergoes a major merger. The last significant 
merger event occurs at z ~ 0.4 (a — 0.71), when the system 
accretes the last 20% of its final mass. After this the system 
remains relatively undisturbed and by z = it is close to 
virial equilibrium. The virial radius, also shown in Figure 3, 
changes by less than 7% after z ~ 0.4. The mass in the 
inner regions of the halo is assembled much earlier. Half of 
it is already in place by z w 5 (a « 0.17) and after z ~ 1 
substantial fluctuations occur only during major mergers. 
(See the triangles in Figure 3, which track the mass in the 
innermost 20 (physical) kpc.) 

2.5 The Analysis 

We focus our analysis on the spherically averaged mass pro- 
file at z — 0. This is measured by sorting particles in distance 
from the centre and binning them in groups of 100 particles 
each. The cumulative mass within these bins, M(r), is then 
used to co mpute the circular velocity profile of each halo, 
V c (r) = v/ GM(r)/r, and the cumulative density profile, 
p(r) = 3 M(r) / Airr 3 , which we shall use in our analysis. 

It is important to choose carefully the halo centre, es- 
pecially since the halos are not spherically symmetric. The 
centre of each halo is determined using an iterative technique 
in which the centre of mass of particles within a shrinking 
sphere is computed recursively until a convergence criterion 
is met. At each step of the iteration the centre of the sphere 
is reset to the last computed barycenter and the radius of 
the sphere is reduced by 2.5%. The iteration is stopped when 
a specified number of particles (typically either 1000 parti- 
cles or 1% of the particles within the high-resolution region, 
whichever is smaller) is reached within the sphere. Halo cen- 
ters identified with this procedure are quite independent of 
the parameters chosen to initiate the iteration, provided that 
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0.4 0.6 
Expansion Factor 

Figure 3. Evolution of the virial radius of the main progenitor of 
the system, r200i of the mass contained within that radius, Af200j 
and of the mass within the innermost 20 h~ 1 (physical) kpc. Data 
are normalized to values at the present day. For the virial radius 
the ratio of the values in comoving units is shown. The system 
undergoes its last major merger at z ~ 0.66, accretes little mass 
afterwards and is close to virial equilibrium at z = 0. The mass 
within the inner 20 h" 1 kpc is assembled earlier than the rest, 
and is only affected seriously during major mergers. 



the initial sphere is large enough to encompass a large frac- 
tion of the system. In a multi-component system, such as 
a dark halo with substructure, this procedure isolates the 
densest region within the largest subcomponent. In more 
regular systems, the centre so obtained is in good agreement 
with centers obtained by weighing the centre of mass by the 
local density or gravitational potential of each particle. We 
have explicitly checked that none of the results presented 
here are biased by our particular choice of centering proce- 
dure. 



3 THE RELATIONSHIP BETWEEN PARTICLE 
NUMBER, SOFTENING, AND TIMESTEP 

The main goal of this study is to identify the conditions 
under which the structure of simulated halos, in particular 
their circular velocity profile, is independent of numerical 
parameters. We start with a brief discussion of the relation- 
ship between three of the main parameters: the number of 
particles, N, the gravitational softening, e, and the timestep, 
At (§ 3.1). We proceed then (§ 3.2) to verify numerically the 
scalings expected between these quantities through a series 
of runs where the timestep for all particles is kept fixed and 
constant throughout the evolution. 



3.1 Analytic Estimates 

Modeling the formation of dark matter halos with N-body 
simulations entails a number of compromises dictated by 
limited computing resources. The choice of particle num- 
ber, timestep, and gravitational softening may all affect, in 
principle, the reliability of the structure of simulated halos. 
We explore here the various limitations imposed by these 
numerical parameters. The analysis assumes, for simplicity, 
a steady-state system with circular speed, V c (r); enclosed 
mass, M(r) = rV c (r) 2 /G; enclosed particle number, N(r); 
and orbital timescale, t c irc — 2irr/V c . The specific energy 
of a typical orbit at radius r is E(r) w v 2 w V c (r) 2 . 



3.1.1 N and Collisional Relaxation 

When a finite number of particles is used to represent a 
system, individual particle accelerations will inevitably de- 
viate from the mean field value when particles pass close 
to each other. Even when orbits are integrated with perfect 
accuracy, these 'collisions' lead to changes of order unity 
in energy on the relaxation timescale (see, e.g., Binney & 
Tremaine 1987), 



(3) 



trelax ^ N (r) 

tcirc In (r/ e) ' 

Thus energy changes due to two-body effects after integra- 
tion time to are given by 



5E_ 

E 



^rcla 



1/2 



tp ln(r/e) 
tcirc (r) N(r) 



1/2 



(4) 



Two-body effects first become important in the inner core 
of the system. Suppressing these effects is primarily a con- 
dition on the number of particles and depends only weakly 
on e. The timestep, of course, does not appear explicitly in 
this criterion. We shall return to the limitations imposed by 
collisional relaxation in S 5.4. 



3.1.2 Timestep and Integration Accuracy 

Accurate integration of the equations of motion of dark mat- 
ter particles requires a careful choice of the timestep adopted 
to evolve the system. A second-order accurate integration 
with timestep At induces a relative error in position, veloc- 
ity, and energy which scales as 



5r 5v SE fvAty ( At 
— oc — oc — oc I oc ( 

r v E \ r J V tcirc 
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Note that this error depends only on the size of At, and 
that it is independent of N and of e, consistent with our 
assumption of a smooth, collisionless system. 

If errors on subsequent timesteps add incoherently, then 
the error at the end of a total integration time, to, is 



6E ( t 
~~E K V At 



V2 / At 
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(6) 



For a given At, then, we expect orbits to be reliably modeled 
only at radii exceeding a certain value r™ defined by, 
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3.1.3 Timestep and Discreteness Effects 

Finite- TV systems are not smooth, and errors in the integra- 
tion will also occur during close encounters between par- 
ticles. The effects of such encounters will be incorrectly 
treated by the simple integrators used in PKDGRAV and 
GADGET whenever the predicted separation at mid-step be- 
tween a particle and a near neighbor satisfies |s| = s < vAt. 
The error in velocity at the end of the step induced by this 
'unexpected' encounter is, then, 



PKDGRAV: At = constant: N ste = 1 28 3 : c=1.25 kpc/h 



8v ~ At 



Gms 



(s 2 + e 2 ) 3 / 2 



(8) 



assuming Plummer softening. Such encounters occur with 
probability 



2 P 

p(s) ds ~ Airs ds — 
m 



<Is 



Gm/i 



(9) 



where p is the mean matter density at the point of encounter, 
and m is the particle mass. The maximum possible size of 
this error is 

GmAt 



(Sv) i 



(10) 



The average velocity change obtained by integrating eq. 8 
over the particle distribution is just that due to the mean 
density field of the system. However, averaging the specific 
energy change over the discrete particle distribution gives a 
positive second-order contribution in excess of that expected 
along the mean-field orbit. For a single step, 



5E. 



GmAt 2 



(11) 



where the simplification arises because the integral is dom- 
inated strongly by contributions at s ~ e. After integration 
time to the total energy change is then, 
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For a given At, then, we expect orbits to be reliably modeled 
at radii larger than a certain r CO nv defined by the following 
condition, 



to 



(At\ 1/2 (Gm/e) 1/2 
V to J V c {r conv ). 



(13) 



Since V c does not change dramatically with radius in CDM 
halos, we see by comparing eq. 7 with eq. 13 that, in the 
presence of discreteness effects, the number of timesteps re- 
quired for convergence increases as e -1 . Economy reasons 
thus dictate the use of large softenings to minimize the 
number of timesteps. On the other hand, large softenings 
compromise the spatial resolution of the simulations. These 
competing effects suggest the existence of an 'optimal' soft- 
ening choice, e op t, which maximizes resolution whilst at the 
same time avoiding discreteness effects and thus minimizing 
the number of timesteps required. We turn our attention to 
the softening next. 




Figure 4. Circular orbit timescale as a function of radius for a 
series of runs with constant timestep. All runs have 128 3 particles 
within the high-resolution box, e = 1.25 h —1 kpc (shown with 
a dotted vertical line), and have been run with PKDGRAV. The 
total number of timesteps used in each run increases from the 
top down, from A?At = 100 to N/\t = 6400 for the dashed curve 
at the bottom. From top to bottom, arrows mark the smallest 
radius where convergence, relative to the smallest-timestep run, 
is achieved in each case. 



individual particle is roughly a me = Gm/e 2 , where m de- 
notes the particle mass. It is useful to compare this with the 
minimum mean field acceleration, which occurs at the outer 
(virial) radius of the system, a m i n ~ GAi^oo /Hmo • The con- 
dition a me < a m i n sets a lower limit to the softening needed 
to prevent strong discreteness effects, 



?"200 
V -/V20( 



(14) 



where JV200 = M2oo/m is the total number of particles 
within r2oo- When this condition is satisfied, discreteness 
causes only small changes in particle accelerations, and so 
does not significantly affect the timestepping in integration 
schemes with an acceleration-based timestep criterion. 

Note that this condition is typically more restrictive 
than the usual requirement that large-angle deflections be 
prevented during two-body encounters. The latter is given 
by t > £26 = Gm/a 2 , where a is the characteristic ve- 
locity dispersion of the system (White 1979). Since a 2 ~ 
GAf20o/2r2oo = Gm A?2oo/2 r2oo, then this condition re- 
quires that forces be softened on scales smaller than 62b ~ 
2r2oo/^V20o, which is usually smaller than e a cc- We shall de- 
termine the relationship between e acc and the 'optimal soft- 
ening' e op t referred to in § 3.1.3 empirically in § 3.2. 



3.1.4 Softening and Discreteness Effects 

When accelerations are softened, the maximum stochastic 
acceleration that can be caused by close approach to an 



3.2 Runs with Constant Timestep 

In order to validate the scalings derived in the previous sub- 
section and to determine empirically the optimum values 
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of the softening and timestep we have carried out a series 
of convergence tests where the timestep has been kept con- 
stant and is shared by all particles. Disabling the multi- 
timestepping capabilities of the codes allows us to concen- 
trate on the role of the timestep size, rather than on the 
virtues or shortcomings of scaling it in various ways for in- 
dividual particles. 

The structure of the dark matter halo chosen for our 
study at z = is illustrated in Figure 4, where we show the 
circular orbit timescale, t c i rc (r) = 2nr/V c (r), as a function 
of radius. Timescales are measured in units of £ c irc(r20o) = 
0.2-kHq 1 , which is of the order of the age of the universe, 
to- Radii are measured in units of the virial radius, r2oo ~ 
205 h~ kpc. The gravitational softening, shown by a vertical 
dotted line, was kept constant in these runs,which had 128 3 
high-resolution particles (~4x 10 5 within r2oo) and were 
run with PKDGRAV. The innermost point plotted in each curve 
corresponds to the radius that contains 100 particles. From 
top to bottom, the curves in Figure 4 illustrate how the mass 
profile of the simulated halo changes as the total number of 
timesteps increases, by successive factors of 2, from N&t = 
100 (top curve) to TVa* = 6400 (bottom dashed curve). 

The halo becomes more centrally concentrated as NAt 
increases, and approaches a 'converged' structure for NAt > 
3200. Runs with fewer timesteps than this still converge to 
the right mass profile but at increasingly larger radii. It is 
interesting to explore how the radius where convergence is 
achieved, r conv , depends on the number of timesteps. We find 
r C onv by identifying the radius at which systematic depar- 
tures greater than 10% in the circular timescale profile first 
become noticeable, gauged against the run with the largest 
number of timesteps. This is easily read off the profiles pre- 
sented in Figure 4. Arrows in this figure indicate t c i IC (r conv ) 
for each choice of NAt ■ 

Filled circles in Figure 5 show the converged timescales 
thus determined as a function of the size of the timestep. 
Converged circular times follow closely the Ai 5//6 depen- 
dence expected from eq. 7, suggesting that the choice of 
softening in this series is such that discreteness effects are 
negligible. This is perhaps not surprising, as e = 1.25 kpc/h 
is about 4 times larger than the lower limit estimated in 
eq. 14, e acc = r 20 o/VN^ = 0.32 kpc/h (for iV sbox = 128 3 ). 

Solid squares in Figure 5 (left panel) correspond to the 
same exercise carried out for several choices of the soften- 
ing when the number of high-resolution particles is reduced 
to 64 3 . For the same softening, e = 1.25 kpc/h, achieving 
convergence with 64 3 particles requires significantly smaller 
timesteps than with 128 3 , as expected since discreteness ef- 
fects become more important as the number of particles is 
reduced. It is also clear from Figure 5 that the dependence 
of tcirc(T'conv) on At changes as e decreases, shifting gradu- 
ally from Ai 5//6 to At 1 / 2 . This transition is precisely what 
is expected from the analytic estimates in § 3.1 (see eqs. 7 
and 13). 

There is further supporting evidence in Figure 5 for 
the validity of the analytic estimates. Consider for exam- 
ple the right panel, where we present the results of runs 
with 32 3 high-resolution particles. The trends are similar to 
those in the left panel, but the transition to the discreteness- 
dominated regime (t c irc(r C onv) x At 1 / 2 ) occurs for even 
larger values of e. 

It is possible to use these results to estimate the soften- 



ing above which discreteness effects become unimportant for 
the various series. From Figure 5, we find that, for 64 3 parti- 
cles, this 'optimal' softening is somewhere between 2.5 and 5 
kpc/h, while for 32 3 particles it is of order ~ 10 kpc/h. Our 
128 3 runs suggest that e op t ~ 1.25 /i -1 kpc for this series. 
The optimal softening appears thus to scale with N just as 
suggested by our discussion of eq. 14. The simple empirical 
rule, 

~ a 4r 2 oo M c s 

£opt " 4eacc = 7^w' (15) 

appears to describe the numerical results well. 

The reason why e opt is about a factor of 4 larger than 
e acc is likely related to the fact that, when softening is 
chosen to optimize results for halos at z = 0, the choice 
is not optimal for their progenitors at earlier times. In- 
deed, r 200 (M,z) oc (n(z)/n ) 1/:i M 1/3 (l + z)-\ which im- 
plies that, for softenings fixed in comoving coordinates, 
e/e op t(N, z) oc N(z) 1 / 6 . Small-iV progenitors thus have 
smaller softenings than optimal and may be subject to dis- 
creteness effects. The dependence on the number of particles 
is weak, however, and it is possible that the factor of 4 in 
eq. 15 may act as a 'safety factor' to ensure that discreteness 
effects are negligible at all times. 

A number of other predictions from the analytic scalings 
presented in § 3.1 are also confirmed by the data in Figure 5. 
For example, when discreteness effects dominate, converged 
timescales are expected to scale as e -1 ^ 2 (eq. 13). This is in 
good agreement with the results of the 32 3 -particle runs; for 
given timestep, t z i IC (r conv ) is seen to increase by roughly a 
factor of 2 when e decreases by a factor of 4, from e = 2.5 
to 0.625 kpc/h. 

Finally, the analytic estimates suggest that the timestep 
choice should be independent of N and e when discreteness 
effects are unimportant. This is also reproduced in the sim- 
ulation series: for e ;> e op t, all runs, independent of N, lie 
along the same dotted line that delineates the 

/ At\ 5/6 

tcirc(r C onv) ~ 15 ^ — J tcirc ( J"200 ) (16) 

scaling. This confirms that the size of the timestep is the 
most important variable when discreteness effects are unim- 
portant; roughly 400, 7000, and 110000 timesteps are needed 
to resolve regions where i c i rc is, respectively, w 10%, 1% and 
0.1% of the orbital timescale at the virial radius. 

3.3 Convergence and integrator schemes 

So far these conclusions are based on runs carried out with 
PKDGRAV. Are they general or do they depend on the partic- 
ular choice of integrator scheme? We have explored this by 
performing a similar series of constant-timestep runs with 
GADGET, which uses a different integrator (§ 2.1 ). There is 
another difference between the GADGET series and the one 
carried out with PKDGRAV: GADGET integrates the equations 
of motion using the expansion factor, a, as the time variable. 
Constant-timestep runs carried out with GADGET were there- 
fore evolved using a fixed expansion-factor step, Aa. Com- 
paring GADGET and PKDGRAV runs with the same total number 
of steps, GADGET takes shorter time steps than PKDGRAV at 
high-redshift, longer ones at moderate z and similar ones at 
z^0. 
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Figure 5. Circular orbit timescale at the smallest 'converged' radius (as illustrated by arrows in Figure 4) as a function of the timcstcp 
for PKDGRAV runs. Left panel shows the results for A^ox = 64 3 , right panel for -/V s t, ox = 32 3 . In both panels the results corresponding to 
A^sbox = 128 3 are shown with filled circles. As the timestep decreases, radii where the orbital timescales are shorter become well resolved. 
This 'saturates' when the radius becomes comparable to the softening and flattens the curves horizontally in some cases. Note also that as 
the softening is reduced the number of timesteps required for convergence increases substantially. Refer to text for a thorough discussion. 



We compare the results of the two series in Figure 6, 
where we plot, at z — 0, the radii containing various mass 
fractions of the halo as a function of the number of timesteps, 
A/a«- The three series shown correspond to runs with 128 3 
high-resolution particles; two were run with GADGET and one 
with PKDGRAV. The choice of softening in each case is indi- 
cated in the figure labels. The four radii shown contain, from 
bottom to top, 0.025%, 0.2%, 1.6%, and 12.8% of the mass 
within r2oo, respectively. 

Convergence is approached gradually and monotoni- 
cally in PKDGRAV runs (solid circles in Figures 5 and 6). For 
A<At ~ 3200 convergence is achieved at all radii containing 
more than ~ 100 particles; fewer timesteps are needed to 
converge at larger radii, as discussed in the previous subsec- 
tion. 

Convergence also occurs gradually, but not monotoni- 
cally, in the case of GADGET (solid squares and triangles in 
Figure 6). For the same number of steps, GADGET results typ- 
ically in mass profiles that, near the centre, are more con- 
centrated than PKDGRAV's, as shown by the systematically 
smaller radii that contain the same mass fraction. The effect 
is particularly noticeable for A^a* ~ 800, when the central 
density profile is actually steeper than the 'converged' result 
achieved for N A t ^ 3200. 

Further runs with different softenings and numbers of 
particles suggest that the presence of these 'cuspy cores' 
in systems evolved with poor time resolution is inherent to 
Aa ^constant GADGET runs, and not just a fluke. On the 
other hand, the artificial cusps only occur in regions well 
inside the convergence criterion derived from the PKDGRAV 
series. Outside the convergence radius delineated by eq. 16 
both GADGET and PKDGRAV results appear safe: one may con- 



clude that GADGET and PKDGRAV require approximately the 
same number of timesteps to resolve the whole system. 

To summarize, the central densities of systems evolved 
with poor time resolution may be over- or under-estimated, 
depending on the integrator scheme adopted. Such sensitiv- 
ity to the integrator scheme emphasizes the vulnerability of 
the central regions to numerical artifact and the need for 
detailed convergence studies such as the one presented here 
before firm conclusions can be reached regarding the inner 
density profiles of CDM halos. 

3.4 Summary 

The agreement presented above between numerical results 
and analytic estimates gives us confidence that it is possible 
to achieve convergence in the mass profiles of simulated dark 
halos down to scales which contain as few as 100 particles or 
where the gravitational softening starts to dominate. A few 
prescriptions for an efficient and accurate integration seem 
clear: 

• choose gravitational softenings so that e > e opt = 
4 rwo / VN200 (eq. 15) to minimize the number of timesteps 
needed, and 

• regard as converged only regions where circular orbit 
timescales exceed w 15 (At/to) s ^ 6 £circ(r20o) (eq. 16). 

One problem with these prescriptions is that, in a large 
cosmological N-body simulation, where systems of different 
mass and size form simultaneously, it is possible to choose 
optimal values of the numerical parameters only for systems 
of roughly the same mass. Also, resolving the inner density 
profiles, where orbital timescales can reach a small fraction 
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Figure 6. Radii enclosing various numbers of particles as a func- 
tion of the total number of timesteps. Results shown correspond 
to runs with 128 3 particles in the high-resolution box and were 
run with PKDGRAV and GADGET, as labeled. PKDGRAV runs approach 
convergence progressively and monotonically. On the other hand 
with poor time resolution GADGET can produce artificially dense 
'cuspy' cores, most noticeable when TV^t = 800. These 'cuspy 
cores' seem to be inherent to the integrator and stepping schemes 
chosen in GADGET. Note, however, that both codes need approxi- 
mately the same number of timesteps for full convergence. 



of the age of the universe, may prove impractical with a 
constant timestep, as the number of timesteps is then dic- 
tated by the densest region of the system, which may contain 
only a small fraction of the total number of particles. It is 
therefore important to learn how the structure of simulated 
dark halos is affected when non-optimal choices of numeri- 
cal parameters are made as well as when multi-timestepping 
integration techniques are adopted. We turn our attention 
to these topics in the following sections. 



4 ADAPTIVE MULTI-STEPPING 
TECHNIQUES 

In order to improve efficiency, many cosmological N-body 
codes use individual timesteps that can vary with time and 
from particle to particle. This allows the time integration 
scheme to adapt spatially so as to achieve high accuracy 
across the whole body of non- linear structures. The two 
codes used in this study, PKDGRAV and GADGET, can use indi- 
vidual timesteps, although, as discussed in § 2.1, they differ 
significantly in the choice of integration scheme. 

Evaluating the efficiency gain is not straightforward, 
since computing resources in most parallel environments do 
not scale in simple ways with the total number of particles 
and of timesteps, and the latter is ill-defined when individ- 
ual adaptive timesteps are adopted. We shall assume, for 
simplicity, that the bulk of the computational work is in- 



vested in computing individual accelerations ('forces'), and 
shall deem efficient timestepping choices that achieve 'full 
convergence'^ whilst minimizing the total number of force 
computations, JVftot- 

For the integrators used in PKDGRAV and GADGET forces 
are computed once every time the position (or velocity) of a 
particle is advanced, so that Aa* = Af to t /N can be thought 
of as the average number of timesteps in a run. Af tot is an 
imperfect measure of the total computational work, since 
it neglects the overhead that stems from tree construction, 
neighbor searching (if required by the timestepping choice) , 
synchronization, and communication between nodes, but is 
nonetheless a useful guide for assessing the efficiency of var- 
ious timestepping techniques. 



4.1 Comparison of timestep criteria 

GADGET allows for five different ways of setting the timestep, 
and we have explored extensively four of them. Our main re- 
sults are illustrated in Figure 7, which is analogous to Fig- 
ure 6 but for runs with 32 3 high-resolution particles. The 
radii shown enclose 1.6%, 3.2%, 6.5%, 12.9%, and 25.8% 
of the mass within the virial radius, respectively, and are 
shown as a function of the timestep parameter, r\ (§ 2.1.1). 
We adopt for this series a softening of 7 h~ x kpc, close to the 
'optimal' value for this number of particles (see Table 1). For 
convenience, we have scaled r\ by an arbitrary factor / (listed 
in the labels of Figure 7) chosen so that, for given / 77, all 
runs in this figure incur approximately the same total num- 
ber of force computations. CPU consumption is lowest for 
EpsAcc and VelAcc, ~ 25% higher for SgAcc, and highest (by 
~ 60%) for RhoSgAcc because the neighbor search required 
by the latter two criteria imposes a significant overhead. 

The main conclusion to be drawn from Figure 7 is that 
all timestepping choices appear to converge for approxi- 
mately the same value of / 77 <; 0.2 or, equivalently, for the 
same AW For fr) ^ 0.2, N ftot « 2.2 x 10 7 , which implies 
that on average a minimum of ~ 650 timesteps is required for 
full convergence. This is comparable to the number of con- 
stant timesteps needed for full convergence (see Figure 5). 

For / 77 > 0.2, deviations from convergence are obvious 
in all cases. Deviations are monotonic in the case of PKDGRAV 
and RhoSgAcc runs; densities at all radii increase gradually 
as the timestep decreases and converge for / 77 < 0.2. On the 
other hand, the behaviour of the inner mass profile in the 
case of other criteria is clearly non-monotonic: the central 
shells dip well below the converged value before bouncing 
back to convergence as / 77 approaches 0.2. This is reminis- 
cent of the artificially cuspy cores discussed in § 3.3, but it 
seems to affect radii well beyond the softening. 

Note that these artificial 'cuspy cores' affect runs with 
GADGET'S EpsAcc criterion as well, which is formally the same 
as used in PKDGRAV. The monotonic approach to convergence 
seen in PKDGRAV runs thus suggest that the presence of 'cuspy 
cores' in runs with poor time resolution is an artifact related 
to GADGET'S integrator scheme rather than to the timestep- 
ping choice. Artificially cuspy cores are an undesirable fea- 



^ We use the term 'full convergence' when it extends down to 
the scale containing as few as 100 particles or the gravitational 
softening, whichever is larger. 
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Table 1. Properties of the simulated halo. 



Label r 2 oo V200 M 2 oo W sbox 7V 2 oo e acc e opt 

[h- 1 kpc] [km s- 1 ] [10 10 M Q ] [h- 1 kpc] [ft" 1 kpc] 

Halo 1 205 205 200 256 3 3.17 x 10 6 0.12 0.46 

128 3 3.97 x 10 5 0.33 1.30 

64 3 4.96 x 10 4 0.92 3.68 

32 3 6.20 x 10 3 2.60 10.4 
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Figure 7. As Figure 6, but for radii as a function of the time- 
stepping parameter, r\. The number of timesteps decreases linearly 
with rj (§ 2). The values of n shown in the Figure have been scaled 
by an arbitrary factor, /, so that for given / rj all runs have similar 
number of total force computations, iVf tot . Values of / are given 
in the figure labels. 

ture in large cosmological simulations, because dense cores 
may survive the hierarchical assembly of structure and lead 
to artifacts in the density profiles of systems formed by the 
merger of affected progenitors. This kind of subtle artifact 
again demonstrates that careful convergence studies of the 
kind presented here are needed to guarantee that the inner 
mass profiles of dark matter halos can be robustly measured 
in N-body simulations. 

4.2 The Dependence on Softening 

According to the analysis presented in § 3.1, the timestep re- 
quired for convergence is independent of the softening when 
discreteness effects are unimportant (i.e., when e J> e op t, 
see eq. 7) but should become increasingly short as e de- 
creases below the optimal value (see eq. 13). Since optimal 
softenings can only be adopted for systems of roughly the 
same mass in a large cosmological simulation, optimizing the 
choice for massive clumps leads to less-than-optimal soften- 
ings in low-mass halos. For such systems, keeping the val- 



ues of rj found to give convergence in the last subsection 
(i.e. f rj ~ 0.2 with the values of / given in Figure 7) may 
not guarantee convergence unless the timestepping criterion 
scales appropriately with softening. For fixed rj, timesteps 
decrease as e 1//2 in PKDGRAV and for the EpsAcc criterion 
of GADGET, but are unchanged as the softening decreases 
for the other GADGET criteria. 

The effects of this are illustrated in Figure 8, which 
shows the result of adopting f rj ~ 0.2 whilst gradually re- 
ducing the softening to values almost two orders of magni- 
tude below optimal. For RhoSgAcJ (solid triangles), f rj = 
0.15 seems appropriate for e close to or slightly smaller than 
fopt ~ 10 h^ 1 kpc, but an artificially low density core clearly 
develops for softenings well below the optimal value. This 
behaviour is not seen in the case of EpsAcc, where conver- 
gence appears firm even for values of e approaching the large 
angle- deflect ion limit, e 2 t- 

We emphasize that this does not signal a failure of the 
RhoSgAcc criterion; rather, it implies that the timesteps cho- 
sen with f rj ~ 0.2 are not short enough to achieve con- 
vergence when e <C e opt . Indeed, choosing RhoSgAcc and 
f rj ~ 0.2(e/e O p t ) 1//2 for small softenings eliminates the arti- 
ficially low density core shown in Figure 8 at a cost in total 
number of timesteps, Nf to t/N, not very different from that 
required by EpsAcc. This demonstrates clearly the need to 
take smaller timesteps when e <C e opt . 

How should timesteps scale with e? Eq. 13 suggests 
a linear dependence when discreteness effects dominate, 
At oc e, although the firm convergence seen for EpsAcc in 
Figure 8 indicate that a gentler dependence, At oc e 1 ^ 2 , 
may actually suffice. This is because the actual individ- 
ual timesteps in this criterion are determined by the ratio, 
(e/cii) 1 / 2 , and accelerations are high during close encoun- 
ters when softenings are small. As a result, the 'effective' 
size of EpsAcc timesteps scales roughly linearly with soften- 
ing when e <C e op t. We have verified this by comparing the 
'maximum' number of timesteps, defined by the total num- 
ber of timesteps taken by a hypothetical particle which, at 
all times, has the minimum timestep of all particles in the 
system, with the minimum number of constant timesteps re- 
quired for convergence (see § 3.2). The agreement is quite 
good. 



II For simplicity, we discuss here only RhoSgAcc; similar results 
apply to VelAcc. 
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Figure 8. Radii enclosing various mass fractions measured at 
z = in our 32 3 simulations as a function of the gravitational 
softening scale length, e. Pairwise interactions become Newtonian 
at distances exceeding 2e. The virial radius of the halo is T200 = 
205 h~ 1 kpc and the total number of particles within this radius is 
N(r20o) ~ 6200. Note that for PKDGRAV and EpsAcc runs the mass 
profile is independent of softening for e < 6h~ 1 kpc, provided 
that the softening remains larger than £2i>> the minimum needed 
to prevent large-angle deflections during particle collisions (§ 3.1). 
For RhoSgAcc, the choice f r) = 0.15 leads to convergence for e > 
Ih^ 1 kpc, but results in large deviations for smaller softenings. 
See § 4.2 for details. 

4.3 Adaptive versus Constant Timestep 

Finally, we investigate the computational gain/loss associ- 
ated with adopting a constant or adaptive time stepping 
technique when a criterion such as EpsAcc is selected. Again, 
we shall assume that the bulk of the computational work 
is invested in computing individual accelerations, although 
this measure neglects the cost of tree construction. Ordinar- 
ily, tree-making contributes a small fraction of the CPU bud- 
get, but this is not necessarily the case in multiple timestep- 
ping schemes when a full tree structure is recalculated every 
time particles in the smallest time bin are advanced. This is 
the case in the version of PKDGRAV that we tested. GADGET, 
on the other hand, recomputes trees only after a certain 
number of interactions have been computed (§ 2.1), so the 
comparison is not straightforward. 

We have chosen for the comparison maximally- 
converged PKDGRAV runs, i.e., those requiring the minimum 
number of timesteps for full convergence. The main con- 
clusion may be gleaned from Table 2, where we list the 
total number of force computations, Nf tot , for runs with 
32 3 high-resolution particles** and three different choices 

For ease of comparison, we have not reduced in this series the 
number of high-resolution particles through the 'amoeba' proce- 
dure described in the Appendix for runs listed in this Table 2. 



for the gravitational softening; e = 10/i _1 kpc (~ e op t), as 
well as e = 2.5 and 0.625 /i -1 kpc. The number of constant 
timesteps needed for full convergence depends sensitively on 
softening, as discussed in § 3; N&t climbs from 800 to 25600 
as e decreases from 10 to 0.625 h~ kpc. The total number 
of force calculations is directly proportional to NAt, and in- 
creases from 2.6 x 10 7 to 8.4 x 10 8 . 

Table 2 shows that, when adaptive multiple timesteps 
are allowed, the total number of force calculations needed 
is comparable when e ~ e op t, but far fewer when the soft- 
ening is well below the 'optimal' value. This demonstrates 
that the small timesteps required when the softening is well 
below the 'optimal' value are only needed briefly by a small 
subset of particles undergoing close encounters. Adaptive 
multi-stepping schemes vastly outperform the fixed timestep 
approach when e <C e opt . 

4.4 Summary 

To summarize, we find that all timestepping criteria we 
have considered can deliver convergence at comparable cost. 
However, the EpsAcc criterion is the one that suffers least 
from overheads related to computing values for individual 
timesteps, and thus appears to be the most efficient of the 
criteria explored in this study. We emphasize, however, that 
this choice is primarily empirical; further investigation may 
very well lead to better and more efficient alternatives than 
any of the ones considered here. 

Further, for softenings close to the 'optimal' value, 
the computational gain that results from adopting multi- 
stepping schemes is rather modest, especially considering 
that the implementation of multi-stepping incurs a non- 
negligible cost in terms of memory usage and bookkeeping. 
Smaller softenings increase the importance of discreteness 
effects and lead to integrations with very small timesteps 
dictated by occasional encounters. Multi-stepping schemes 
are strongly favored under these circumstances. 



5 THE ROLE OF OTHER NUMERICAL 
PARAMETERS 

Proper convergence requires, of course, that appropriate 
choices be made for all relevant parameters. We now turn to 
the analysis of the separate role of other numerical param- 
eters. Unless explicitly stated, we will undertake the anal- 
ysis of each parameter using only runs for which all other 
parameters take 'converged' values. This can only be done 
after a large parameter space search since the effects of com- 
binations of some parameters may be subtle. For example, 
a timestep that is adequate for some gravitational softening 
may be inadequate when the softening is substantially mod- 
ified. Because of this restriction, the results in the following 
subsections contain, for clarity, only a small fraction of all 
runs performed. 

5.1 The Gravitational Softening 

Large cosmological simulations generally use a single par- 
ticle mass and thus resolve systems of different mass with 
different numbers of particles. This implies that it is pos- 
sible to choose 'optimal' values of the softening only for a 
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Table 2. Properties of maximally-converged runs (PKDGRAV). 
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small range of halo masses, since e op t oc V2oo /N 2 q oc N^qq ■ 
This may not be too restrictive for the resimulations we dis- 
cuss here, since they focus on one system at a time, but it 
does affect significantly large cosmological simulations. If, 
for example, an optimal softening choice is made for the 
most massive system expected to form at, say, z = , it will 
be smaller than the optimal value for less massive systems 
present at the same time (see eq. 15). How are their mass 
profiles affected and what regions in such systems may be 
considered converged? 

To address this question, we have undertaken a large 
series of simulations where the softening, e, was varied sys- 
tematically while choosing 'converged' values of all other 
parameters. We have explicitly checked that, for example, 
doubling or halving the number of timesteps (or the ini- 
tial redshift) has no appreciable effect and that, for given 
number of particles, the results discussed in this subsection 
depend only on e. 

We show the results of this series in Figure 8, where 
radii enclosing various mass fractions are shown as a function 
of e in simulations with 32 3 high-resolution particles. Since 
?"200 ~ 205 h~ kpc, the radii shown in Figure 8 probe a 
large fraction of the halo's radial extent, between 4% and 
22% of the virial radius. For this system, e2& = 0.066 h~ kpc 
« 3.2 x 10" 4 r 20 o and e opt ~ 10 h' 1 kpc « 4.9 x 10" 2 r 20 o- 

As Figure 8 shows, the mass profiles obtained with the 
two codes agree to better than 20% (i.e., to better than 
10% in circular velocity) even for radii containing as few as 
100 particles. Full convergence is achieved for a wide range 
of softening scales, provided that £26 < e <, 6ft" 1 kpc. The 
mass profiles are essentially unchanged even as the softening 
is varied by almost two orders of magnitude. 

A second important point to note in Figure 8 is that 
for e ~ 12 ft -1 kpc (only slightly larger than e op t) the profile 
deviates from the converged one even as far out as 60 /i -1 
kpc; i.e., more than 5 times the softening length. This con- 
trasts with the results for e ~ 6 ft" 1 kpc, where the mass 
profile appears to have converged down to almost one soft- 
ening length scale. Clearly, assuming that mass profiles are 
affected out to a certain multiple of the softening length is 
an oversimplification that is not supported by these results. 

What determines the smallest converged radius for a 
given softening length scale? Since softenings introduce a 
characteristic acceleration on small scales, it is instructive to 
consider the mean acceleration that particles experience as a 
function of the distance from the centre of the system. This 
radial acceleration profile, a(r) — GM(r)/r 2 = V c 2 (r)/r, is 
shown in Figure 9 for two series of runs where the grav- 
itational softening has been varied systematically by two 




-3 -2 -1 

L °g r / r 2oo 



Figure 9. Spherically-averaged 'acceleration' profiles (V c 2 (r)/r) 
for 64 3 and 32 3 runs, shown for several choices of the softening 
scalelength, e. The dotted line corresponds to the acceleration 
profile of an NFW model with concentration c = 10. The ver- 
tical arrows denote the value of the softening parameter, e, for 
each run. The profiles line up, from bottom to top, in order of 
decreasing e. As e approaches ~ 0.01 r200, the acceleration pro- 
files converge to a solution similar to the fiducial NFW curve. 
Profiles significantly affected by the softening deviate from the 
converged result at a radius where the acceleration matches the 
characteristic acceleration associated with the circular velocity of 
the halo, V200, and e: a e = Xe^rjo/ 6 ' w h-h Xe ~ 0.5. The values 
of a t corresponding to each adopted value of e are shown by the 
horizontal arrows. 

orders of magnitude. The values of the softening in each run 
are shown with small vertical arrows near the bottom of the 
figure. Solid and dashed curves correspond to runs with 64 3 
and 32 3 particles in the high-resolution box, respectively. 
As the softening is decreased from e ~ 0.1 r2oo by successive 
factors of two, the acceleration profiles become steeper and 
converge to a unique profile for e < 0.03 r2oo ~ 6 ft -1 kpc, 
as shown in Figure 9. The convergent profile is well approx- 
imated by an NFW profile with c = 10, shown by a dotted 
line in Figure 9. 

We note two interesting features of the acceleration pro- 
files shown in Figure 9. The first is that the effects of soft- 
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ening on the acceleration profile depend rather weakly on 
the number of particles used; for given e, the profiles corre- 
sponding to runs with 32 3 and 64 3 particles agree reason- 
ably well, and they approach the same 'converged' profile for 
e < 0.03 r 2oo- The second feature is that acceleration pro- 
files deviate from the 'converged' profile near the centre for 
larger values of the softening. Interestingly, deviations oc- 
cur at radii where the acceleration exceeds a 'characteristic' 
acceleration, 



E=1.25 kpc/'h: EpsAec 



(17) 



which depends only on the circular velocity of the halo and 
on the value of the softening adopted. This characteristic 
acceleration is shown (for Xe ~ 0.5) with horizontal arrows 
in Figure 9. The mass profile of a simulated halo becomes 
unreliable for accelerations exceeding a e . 

This result suggests an empirical interpretation of the 
effects of softening on the mass profile of a simulated halo: 
the choice of gravitational softening imposes an effective 
limit on the accelerations that may be adequately repro- 
duced in the system. This is interesting, since for systems 
with density profiles similar to that proposed by NFW, there 
is a maximum acceleration that particles may experience. 
Indeed, a(r) = V c 2 /r tends to a well-defined maximum, 



c 2 /2 



ln(l + c) — c/(l + c) r2oo 
as r approaches zero. If e is such that 



(18) 



(19) 



then it appears to impose no substantial restriction on the 
mass profile. For example, Figure 6 shows that the converged 
mass within ~ l/i -1 kpc appears not to change as e varies 
between 1.25 and 3.5/i _1 kpc. At face value, this would ap- 
pear to imply that the mass profile can be trusted down 
to almost one third of the softening length scale when the 
condition expressed in eq. 19 is satisfied. In order to be con- 
servative, however, we shall hereafter assume that converged 
radii cannot be less than e. 

How does the upper limit on e dictated by this con- 
straint compare with e op t, the minimum needed to prevent 
discreteness effects and minimize the number of timesteps? 
The answer depends on the number of particles, as well as 
on the concentration of the system, and imposes an effective 
lower limit on the number of particles needed to satisfy both 
conditions simultaneously, -/V200 > (2c) 4 /(ln(l + c) — c/(l + 
c)) 2 . Forcw 10, we find that roughly 70, 000 particles within 
the virial radius are needed to carry out a simulation where 
the softening is small enough not to restrict significantly the 
resolution of the inner mass profile and large enough to pre- 
vent discreteness effects from hindering the computational 
efficiency of the calculation. 

To summarize, provided that all other numerical pa- 
rameters are chosen appropriately, the effect of the soften- 
ing on the spherically-averaged mass profile is to impose a 
maximum acceleration scale above which results cannot be 
trusted. The mass profile of a simulated halo converges at 
radii where the mean acceleration does not exceed a char- 
acteristic value imposed by the softening, a(r) = V^(r)/r < 
a e = Xe^ooAi where Xe is empirically found to be ~ 0.5 if 
e is expressed as a spline-softening scalelength. 




Figure 10. As Figure 8, but as a function of the initial rcdshift 
of the simulation. Convergence is seen for Zi > 25 in the 32 3 runs 
and for Zi > 49 in the 64 3 runs. Starting at lower initial rcdshifts 
causes halo mass profiles to develop an artificially low density 
core. 

5.2 The Initial Redshift 

The starting redshift, z;, determines the overall initial ampli- 
tude of density fluctuations in the simulation box. If Zi is too 
low, small scales may already be in the non-linear regime, 
invalidating the assumptions of the procedure outlined in 
§ 2.2. Initial redshifts cannot be chosen to be too high ei- 
ther, since the more uniform the periodic box, the more 
difficult the task of evaluating accurate forces in treecodes 
such as the ones we employ here becomes. A compromise 
must therefore be struck between these competing demands 
and we derive in this section a simple empirical prescription 
that ensures convergence in the mass profiles of simulated 
CDM halos at z = 0. 

Figure 10 shows the radii of various mass fractions (at 
z = 0) as a function of the initial redshift of the simulation. 
Top and bottom panels refer to the same halo, using two 
different particle numbers in the high-resolution box: 32 3 
(bottom), and 64 3 (top). Each curve is labeled by the en- 
closed number of particles. The inner mass profile of the halo 
converges as the initial redshift is increased. Convergence to 
better than 10% at all radii is achieved for 25 < 1 + Zi < 100, 
and even for the highest Zi tested there is no clear departure 
from convergence. We have checked explicitly that this re- 
sult does not depend on the particular time-stepping choice; 
a similar series with the SgAcc criterion gives similar results. 

The data in Figure 10 also suggest that convergence 
may be achieved at lower Zi when 32 3 particles are used 
rather than 64 3 . A possible explanation for this is presented 
in Figure 11, where we plot, for each radial shell, the devi- 
ations from the converged value as a function of the (theo- 
retical) rms mass fluctuation on the smallest resolved mass 
scale at Zi, o(m p , Zi) (m p is the mass of one high-resolution 
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Figure 11. The radii of various mass shells, as in Figure 10, but 
normalized to the 'converged' value of the radius for each shell as 
a function of a(m p , Zi), the linear rms fluctuation on the scale of 
the particle mass at z = Zi. Note that convergence is achieved at 
all radii when tr(m p , zi) < 0.3. 



particle). In terms of this variable, the 32 3 and 64 3 results 
are indistinguishable, showing convergence down to the 100- 
particle mass shell when Zi is chosen so that o~(m p , Zi) <J 0.3. 
This is a simple empirical rule for choosing the starting red- 
shift that we shall adopt hereafter. 

One advantage of this rule is that, for power spectra 
such as CDM, a(m p ) is only weakly dependent on mass on 
small scales, so the initial redshift can be chosen almost 
independently of the number of particles. For example, even 
for the highest number of particles considered in our study 
(iVsbox = 256 3 , m p = 6.5 x 10 5 /i _1 M ) the starting redshift 
condition is satisfied for 1 + zi J> 42, so that 1 + m = 50 
could be safely used for all of our simulations, regardless of 
N. 



5.3 Force Accuracy 

Accurate forces are an obvious requirement for numerical 
convergence, and we investigate here the role of force accu- 
racy parameters in the mass structure of dark halos. This 
is important since treecodes are based on approximate mul- 
tipole expansion-based methods that are vulnerable to in- 
accuracies in the force calculations. Although accuracy can 
always be improved by adopting, for example, stricter node- 
opening criteria, this comes usually at the cost of substantial 
loss in computational efficiency. It is therefore important to 
determine what is the minimum force accuracy needed to 
achieve convergence in order to maximize the efficiency of 
the simulation. 

Force accuracy is controlled in GADGET (in the config- 
uration used in this study, see § 2.1.1) through two main 
parameters: a compile-time flag, -DBMAX, which, if enabled, 




Figure 12. As Figure 8, but for radii as a function of the GADGET 
force accuracy parameter, / acc (ErrTolForceAcc in GADGET'S pa- 
rameter file). Filled squares show results without enabling the 
extra-accuracy flag -DBMAX during compilation. Filled circles show 
results enabling -DBMAX. When this flag is on, the effects of / aC c 
on the mass are mild, and good convergence is achieved even for 
rather large values of / aC c- When -DBMAX is off, / acc < 10~ 3 is 
needed to ensure convergence. 



restricts node opening to a list of cells guaranteed not to con- 
tain the particle under consideration, and by the parameter 
/ acc (named ErrTolForceAcc in GADGET'S parameter file), 
which controls dynamically the updating of the tree-node 
opening criterion (Springel et al. 2001). Figure 12 shows the 
radii of various mass shells in our standard halo as a function 
of / a cc • Filled squares show the results obtained without set- 
ting the -DBMAX option in GADGET. Convergence is achieved 
in this case for quite small values of the accuracy parameter, 
/acc <, 0.003. The reason behind the slow convergence seen 
in Figure 12 appears to be related to rare but substantial 
errors incurred in GADGET'S tree walking procedure when the 
boundaries of open nodes are not guaranteed to exclude the 
particle under consideration (Salmon & Warren 1993). Dis- 
allowing this possibility (i.e., enabling -DBMAX during compi- 
lation) leads to much improved convergence relative to the 
parameter /acc, as can be seen from the filled circles in Fig- 
ure 12. There is almost no systematic trend with / aC c when 
-DBMAX is enabled, even for / acc ~ 1. 

The main effect of enabling -DBMAX is to suppress a tail 
of large errors that, although rare, appear to have a signifi- 
cant effect on the final mass profile. This can be seen in Fig- 
ure 13 where we show the cumulative distribution of errors 
in accelerations computed on a z — snapshot of a simu- 
lation with 32 3 high-resolution particles. Force errors were 
measured by comparing with the result obtained by direct 
summation. Solid and dashed lines give the result of open- 
ing nodes with the 'relative' opening criterion proposed by 
Springel et al. (2001), with and without the -DBMAX option 
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Figure 13. Cumulative error distributions of GADGET'S force com- 
putation for various choices of opening criterion and tolerance pa- 
rameter. We used the particle distribution of the z=0 snapshot of 
a run with 32 3 high-resolution particles and measured force errors 
by comparing to the result obtained by direct summation. Solid 
and dashed lines give the result of opening nodes with the 'rel- 
ative' opening criterion proposed by Springcl ct al. (2001), with 
and without the -DBMAX option (solid and dashed lines, respec- 
tively). Results are shown for tolerance parameters / acc = 0.001, 
0.003, and 0.01 (from left to right). Dotted lines show results for 
the traditional BH-opcning criterion (dotted lines) , with opening 
angles 9 = 0.5, 0.75, 1.0 (from left to right). The inset compares 
the accuracy obtained for all of these choices as a function of 
computational cost. See text for more details. 



(solid and dashed lines, respectively). In each case, results 
are shown for tolerance parameters / acc — 0.001, 0.003, and 
0.01 (from left to right). We also show results for the tra- 
ditional Barnes- Hut opening criterion (dotted lines), with 
opening angles 9 — 0.5 , 0.75, 1.0 (from left to right). 

We have chosen a rather small value of the softening in 
Figure 13 to emphasize graphically the point that a long tail 
of errors may exist when the -DBMAX option is not enabled; 
note, for example, that errors of up to 100% or larger are 
present in this case when / acc = 0.01 (rightmost dashed 
line). Such errors are not present when -DBMAX is on (solid 
lines) . 

The inset compares the accuracy obtained for all of 
these choices as a function of the invested computational 
cost. 'Accuracy' is here taken as the 98% percentile force 
error, and the computational cost is measured in terms of 
the average number of node-particle interactions per force 
evaluation. For a given accuracy, the Barnes-Hut criterion 
results in higher cost than the criterion adopted in GADGET. 

We conclude that enabling -DBMAX and adopting / acc <J 
0.01 is sufficient to study the inner structure of dark matter 
halos. Alternatively, adopting a redshift-dependent Barnes- 
Hut node opening criterion, such as in PKDGRAV, where 9 = 



0.55 is used for z > 2 and 6 — 0.7 for z < 2, seems also to 
give adequate results. 

5.4 The Number of Particles 

The total number of particles is a critical parameter to 
choose when running a cosmological N-body simulation. 
Since the computation time will scale at best linearly with 
N, one must try and use as few particles as possible to 
achieve the goals of the programme. As mentioned in § 1, our 
main goal is to provide robust and accurate measurements 
of the circular velocity (or mass) profile of dark matter ha- 
los down to about the inner 1% of the virial radius. This 
corresponds to ~ 2.2 h^ 1 kpc in the case of the Milky Way 
if its halo has the same circular velocity as the disk. This 
is clearly the minimum resolution required for meaningful 
comparison with observed rotation curves. 

In the preceding discussion we have determined the op- 
timal choice of softening, time stepping, force accuracy, and 
starting redshift required to obtain repeatable and robust 
measurements of the circular velocity profile of a simulated 
CDM halo down to radii containing as few as 100 particles. 
Repeatability and robustness relative to these parameters 
are, of course, necessary conditions for convergence, but we 
must still demonstrate that the results do not depend on the 
total number of particles chosen. 

How many particles must a region contain so that the 
circular velocity (or, equivalently, the mean inner density) 
converges? We use the lessons from the preceding subsec- 
tions to explore the dependence of the mass profile of sim- 
ulated dark halos on the number of particles used. We con- 
sider only runs which meet the requirements discussed pre- 
viously, so that, for each choice of N, we shall only present 
'converged' results relative to other parameters. Our tests 
span an unprecedented range of 512 in particle number, from 
32 3 to 256 3 particles in the high-resolution simulation cubes. 

Our main results are summarized in Figure 14, where 
we show, as a function of the enclosed number of particles, 
the mean inner density contrast measured at various radii 
from the centre of the halo. In this figure, for example, solid 
triangles show the mean inner density contrast measured at 
~ 20% of the virial radius. ^,From left to right, each group 
of filled triangles indicates the results of runs with 32 3 , 64 3 , 
128 3 , and 256 3 particles in the high resolution cube. These 
runs have 6200; 49600; 397000; and 3.2 x 10 6 particles within 
r2oo , respectively. As the number of particles increases fewer 
runs are shown, because of the increasing computational 
cost. At the highest resolution, with 256 3 particles in the 
high-resolution cube, we have completed only one simula- 
tion. This run is comparable to the highest resolution simu- 
lations reported in the literature so far. 

Figure 14 shows a number of important trends. Con- 
sider, for example, the radius corresponding to 2% of r2oo 
(solid circles). In the 32 3 runs, this radius contains 1.6% of 
the halo mass (~ 100 particles). The mass within this radius 
is seen to increase significantly as the number of particles in- 
creases; the density contrast climbs from ~ 1.2 x 10 s (in the 
32 3 runs) to ~ 2.5 x 10 5 (in the 64 3 runs) before stabilizing 
at ~ 3 x 10 5 when iV sbox reaches 128 3 and 256 3 . Clearly, 100 
particles are not enough to trace reliably the mass profile 
of a simulated halo, in disagreement with the conclusions of 
Klypin et al. (2001), who argue that 100-200 particles suffice 
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Figure 14. Mean inner density contrast as a function of the enclosed number of particles in 4 series of simulations varying the number 
of particles in the high-resolution box, from 32 3 to 256 3 . Each symbol type corresponds to a fixed fraction of the virial radius, as shown 
by the labels on the right. The number of particles needed to obtain robust results increases with density contrast, roughly as prescribed 
by the requirement that the collisional relaxation timescale should remain longer than the age of the universe. According to this, robust 
numerical estimates of the mass profile of a halo are only possible to the right of the curve labeled f ro i ax ~ 0.6fo- 



to resolve the inner mass profile when other parameters are 
chosen properly. 

The situation is different for the 1000-particle radius in 
the 32 3 runs, which correspond to about 10% of the virial 
radius. The density contrast within this radius is ~ 2.5 x 10 4 , 
and remains essentially unchanged as the number of parti- 
cles increases by a factor of 512. The data presented in Fig- 
ure 14 thus support the conclusions of Moore et al (1998): 
resolving regions closer to the centre, where the density con- 
trast is higher, demands increasingly large particle numbers. 
Although 300 particles in the 32 3 runs are almost enough to 
resolve the 6% radius, they fall well short of what is needed 
to resolve the much higher overdensities characteristic of the 
1% radius. 

How many particles are needed to resolve a given ra- 
dius? Moore et al. (1998) propose that converged regions 
are delineated by (one-half) the mean inter-particle separa- 
tion within the virial radius, 0.5 (47r/3iV2oo) 7-200, whereas 
Fukushige & Makino (2001) suggest that the innermost re- 
solved radius cannot be smaller than the radius where the 
two-body relaxation time becomes shorter than the age of 
the universe. 

Our results appear to favor the latter interpretation. For 



example, the criterion of Moore et al. would predict that the 
32 3 runs could be trusted down to 4.5% of the virial radius, 
but it is clear from Figure 14 that convergence in this case 
is achieved only for radii beyond 6% of 7-200. On the other 
hand, all simulations can be seen to converge at radii larger 
than the radius where the average collisional relaxation time 
roughly matches the age of the universe. This is shown by 
the (almost vertical) line labeled t rc iax ~ to, where we define 

trdax(r) = N rJVc = V200 N t J_\ ~ 1/2 
£circ(r-20o) 8 In N 7-200/^200 8 IniV \p C vit J 

t c itc(i"20o) ~ to, and N — N(r) is the enclosed number 
of particles. For reference, the curve on the left indicates 
trciax = 0.6 tcirc (7-200) ~ 0.6 to • As shown in Figure 14, the 
density profile converges at radii that enclose enough parti- 
cles SO that trelax(T-) > 0.6 to ■ 

We emphasize that this criterion is mainly empirical, 
and does not necessarily imply that particles in regions 
where the relaxation time is shorter than ~ 0.6 to actu- 
ally evacuate the central regions as a result of two-body 
encounters. Indeed, one would expect the inner mass profile 
to evolve as a result of collisions on the much longer 'evap- 



© 2002 RAS, MNRAS 000, 1-19 



Inner halo structure I: convergence 19 



oration' timescale, t CV a P ~ 136t rc i ax (Binney & Tremaine 
1987), a proposition that finds support in simulations of the 
evolution of isolated equilibrium iV-body systems (Hayashi 
et al 2002). In addition, the heating rate near the centre is 
likely dominated by the presence of substructure rather than 
by particle-particle collisions, complicating the interpreta- 
tion. Our result is thus reminiscent of the work of Weinberg 
(1998), who emphasizes the difficulty of achieving the col- 
lisionless limit in iV-body systems and the possibility that 
fluctuation noise may lead to relaxation effects important 
on all scales. 

Despite this difficulty, it seems clear from Figure 14 that 
resolving density contrasts exceeding 10 6 requires > 3000 
particles within that radius, or over 3 million particles within 
the virial radius. Providing robust numerical predictions of 
the mass structure of cold dark matter halos on scales that 
can be compared directly with observations of individual 
galaxies is thus a very onerous computational task. 



5.5 Optimal parameters - a worked example 

The many considerations discussed in the previous sections 
make the selection of optimal parameters for any given N- 
body run a delicate and complicated business. It may be 
helpful to go through how one might choose optimal pa- 
rameters for a specific calculation, for example a simula- 
tion like the largest one (iV s box = 256 3 ) we consider in this 
paper. This run has ~ 3 x 10 6 particles within the virial 
radius at z = 0, and is the largest we can easily carry 
out with resources currently available to us. Figure 14 and 
the discussion in § 5.4 suggest that this number of parti- 
cles should be sufficient to get converged results down to 
about r conv — 0.005 r2oo- Equation 15 suggests that a soft- 
ening parameter e — 0.0025 r2oo will be near optimal for 
getting an efficient integration almost unaffected by dis- 
creteness effects. As Figure 9 demonstrates, this softening 
is small enough relative to our target r conv that it should 
not compromise the radial structure. For these parameters, 
equation 16 and Figure 4 then show that a single-timestep 
integrator should be able to converge in about 5000 equal 
steps, although we note that this depends on the detailed 
inner structure of the halo, which is what we are trying 
to measure. In practice, a series of runs where the num- 
ber of particles is gradually increased, is desirable to fine- 
tune the choice of timestep. Alternatively, the discussion of 
§ 4 implies that for our preferred multi-timestep integrator 
(EpsAcc) 77 = 0.15 should be small enough to ensure conver- 
gence. The discussion of § 5.2 shows that it should be safe 
to start the integration at Zi = 49. 



6 THE CIRCULAR VELOCITY PROFILE OF A 
ACDM HALO 

Finally, we use the convergence lessons derived above to an- 
alyze briefly the inner circular velocity profile of the ACDM 
halo considered here. The results of 'converged' runs are 
shown in Figure 15. Each profile is shown only for radii con- 
sidered converged according to the criteria discussed above. 
Plotted this way, all profiles, independent of the number 
of particles, seem to agree to within ~ 10% at all radii. 



Converged results: GADGET & PKDGRAV 
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Figure 15. Circular velocity profiles of 'converged' runs with dif- 
ferent number of high-resolution particles. Profiles are only plot- 
ted for radii where the convergence criteria derived in this paper 
are satisfied. Several curves are shown for the cases of 32 3 , 64 3 , 
and 128 3 particles, corresponding to runs where all other numeri- 
cal parameters take converged values. For clarity, a small selection 
of runs have been chosen; those with softenings indicated by the 
small vertical arrows. The convergent profile that emerges for this 
halo is roughly independent of the number of particles and resem- 
bles closely the model proposed by NFW, with c = 10. For this 
halo, steeply-cusped density profiles are disfavoured. The profiles 
labeled 'Moore et al' and 'NFW have been matched at the peak 
of the circular velocity profile. 

The circular velocity increases from the virial radius in- 
wards, reaches a maximum and then drops gradually to- 
wards the centre, following closely the dotted line that repre- 
sents an NFW profile with concentration c = 10. This value 
of the concentration agrees reasonably well with the results 
of NFW and of Eke, Navarro & Steinmetz (2001), who find 
c w 8-9 for a halo of this mass. Near the centre, the profile is 
seen to deviate significantly from the steeply cusped profile 
approaching a central slope of j3 = 1.5 proposed by Moore 
et al. (1998), and agrees better with shallower central slopes 
such as that of the NFW model. 

We emphasize that there is little evidence for conver- 
gence to a power-law density profile near the centre, and 
that the profile keeps getting shallower down to the inner- 
most point that our procedure deems converged. Can our re- 
sults be used to place meaningful constraints on the asymp- 
totic inner slope? At r m i n ~ l/i -1 kpc, the smallest ra- 
dius resolved in our highest-resolution run (iV a box = 256 3 ), 
both the local and cumulative density profiles are robustly 

determined^: p(r min )/p cri t = 9.4 x 10 5 , and p(r min )/p cri t « 

tt Convergence in the local density actually extends to radii 
smaller than the minimum converged radius for the more strin- 
gent cumulative density. 
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1.6 x 10 6 . These values can be combined with the require- 
ment of mass conservation to place an upper limit to the 
inner asymptotic slope of the density profile, /3 < 3(1 — 
p(r m in)/p('"min)) = 1.2. In other words, there is not enough 
mass within r m i n to support a power-law density profile with 
slope steeper than /3 = 1.2. We note that this conclusion 
depends sensitively on our ability to resolve the innermost 
lh^ 1 kpc. If r m i n were just two or three times larger the 
same exercise would not be able to rule out slopes as steep 
as /3 = 1.5. 

In summary, our results argue strongly against the very 
steep central cusps advocated by Moore et al. (1998, 1999), 
Ghigna et al. (1998, 2000) and Fukushige & Makino (1997, 
2001). We are in the process of augmenting our sample of 
halos in order to firm up this conclusion, so will defer a 
detailed analysis of this issue to a later paper in this series. 



7 CONCLUSIONS 

We have performed a comprehensive series of convergence 
tests designed to study the effect of numerical parameters 
on the structure of simulated cold dark matter halos. Our 
tests explore the influence of the gravitational softening, the 
time-stepping algorithm, the starting redshift, the accuracy 
of force computations, and the number of particles on the 
spherically-averaged mass profile of a galaxy-sized halo in 
the ACDM cosmogony. We derive, for each of these param- 
eters, empirical rules that optimize their choice or, when 
those choices are dictated by computational limitations, we 
offer simple prescriptions to assess the effective convergence 
of the mass profile of a simulated halo. Our main results can 
be summarized as follows: 

(i) Timestep and Discreteness Effects. The number of 
timesteps required to achieve convergence depends primar- 
ily on the orbital timescale of the region to be resolved, 
but may also be sensitive to the number of particles and 
the gravitational softening, unless these parameters are cho- 
sen so that discreteness effects are unimportant. This re- 
quires the gravitational softening to be large enough so that 
the maximum acceleration during two-body encounters does 
not exceed the minimum mean field acceleration in the halo, 
e > e»cc = ^oo/VA^oo- Empirically, we find that e w e op t = 
4 e aC c gives good results. When this condition is satisfied, the 
minimum converged radius, r CO nv, is given by the condition 
that the circular orbit timescale should be long compared 
to the timestep, t c irc(j"conv) ~ 15 (At/t ) 5/e icirc^oo)- Sub- 
stantially smaller timesteps are needed if e < e op t- Dark mat- 
ter densities at r < r con v may be under- or over-estimated, 
depending on the integrator and timestepping schemes used. 
For example, constant-timestep GADGET runs develop arti- 
ficially dense, 'cuspy' cores in poorly resolved regions, in- 
dicating that the approach to convergence is not always 
monotonia This emphasizes the importance of comprehen- 
sive convergence tests such as the ones presented here to 
validate the results of numerical studies of the inner struc- 
ture of CDM halos. 

(ii) Fixed Timestep versus Adaptive Multi- Stepping. Of 
the several adaptive, multiple time-stepping criteria that we 
considered, we have found best results when timesteps are 
chosen to depend explicitly on the gravit ational softening 
and on the acceleration, AU = r) ac yj u/a,i, with rjae ~ 0.2. 



Experiments with time-stepping choices that do not include 
explicitly the gravitational softening require the value of the 
corresponding n to be reduced as e is reduced below the 
optimal value in order to obtain convergence. In terms of 
computational cost, we find that multi-time-stepping crite- 
ria significantly outperform the use of a single timestep for 
all particles only for softenings well below the optimal value. 

(iii) Gravitational Softening. The choice of gravitational 
softening is found to impose a maximum acceleration scale 
above which simulation results cannot be trusted. This ac- 
celeration scale appears to depend mainly on the circular 
velocity of the halo and on the gravitational softening scale, 
and is given by a e = \ e V^ooAi with x« ~ 0.5. For given par- 
ticle number, convergence to better than 10% in the mass 
profile is obtained at radii greater than e that also contain 
more than 100 particles and where the acceleration criterion 
is satisfied: a(r) — V c (r) 2 /r <J a e . 

(iv) Starting Redshift. The mass profiles of simulated dark 
halos converge provided that the initial redshift is chosen so 
that the theoretical (linear) rms fluctuations on the small- 
est resolved mass scale, m p (the mass of one high-resolution 
particle) is a(m p ,Zi) < 0.3. Since o(m p ) is a weak function 
of mass on subgalactic mass scales for CDM-like power spec- 
tra, this criterion indicates that a modest starting redshift, 
such as 1 + Zi ~ 50 is appropriate for particle masses as low 
as m p ~ 10 5 ft -1 Mq in the ACDM cosmogony. 

(v) Force Accuracy. The mass profiles of simulated CDM 
halos are quite sensitive to the accuracy of the force calcu- 
lations, and convergence requires care in the choice of node 
opening criteria in the treecodes used in our study. Poor 
force accuracy leads to the development of artificially low 
density cores. In the case of GADGET, for example, we find 
that even occasional large errors in the forces may lead to 
noticeable deviations from converged profiles. To avoid this, 
it is necessary to choose tree-walking parameters that cur- 
tail drastically the tail of the most deviant force calculations, 
however rare. In GADGET this can be achieved by activating 
the compiler option -DBMAX. Using up to hexadecapole terms 
in the node potential expansion and setting a redshift de- 
pendent tree-node opening criterion, as in PKDGRAV, where 
6 — 0.55 is chosen for z > 2 and 8 = 0.7 for z < 2, seems 
also to work well. 

(vi) Particle Number. In order to achieve convergence in 
the mass profile, enough particles must be enclosed so that 
the average two-body relaxation timescale within the region 
is comparable or longer than the age of the universe. We find 
empirically that the condition, t re iax(f) > 0.6 to, describes 
converged regions well. Since t relax IS roughly proportional 
to the enclosed number of particles times the local dynami- 
cal timescale, resolving regions near the centre, where den- 
sity contrasts are high and dynamical timescales are short, 
requires substantially more particles than resolving regions 
more distant from the centre. Of order 3000 enclosed parti- 
cles are needed to resolve regions where the density contrast 
reaches 10 6 . On the other hand, density contrasts of order 
10 4 ' 5 require only 100 enclosed particles for numerical con- 
vergence. Resolving radii of order 0.5% of the virial radius 
in the first case requires of order 3 x 10 6 particles within the 
virial radius. 

For most simulations, the most stringent convergence 
criterion is the relaxation timescale condition on the num- 
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ber of particles. This implies that there is little choice but 
to strive for the largest possible N when studying the inner 
regions of dark matter halos. This limit is dictated by the 
available computer resources. Choosing the optimal soften- 
ing for the adopted number of particles then minimizes the 
number of timesteps needed to achieve convergence down 
to the radius where t Te i a *_(r) > 0.6 to- The precise number of 
timesteps cannot be determined ahead of time, since irciax(^) 
depends on the detailed structure of the halo, which is what 
we are trying to measure. This implies that a series of simu- 
lations where the number of particles is increased gradually 
is advisable in order to ensure that optimal parameters are 
chosen for the highest-resolution run intended. 

We have applied our convergence criteria to a ~ 205 
km s _1 ACDM halo in order to investigate the behaviour 
of the inner slope of the density profile. We find that 
the slope of the spherically-averaged density profile, f3 = 
— dlog(p)/dlog(r), becomes increasingly shallow inwards, 
with little sign of approach to an asymptotic value. At the 
smallest radius that we consider resolved in our highest- 
resolution (256 3 ) simulation (r m i n ~ 1 h^ 1 kpc « 0.005 r2oo), 
the local and cumulative density contrasts are robustly de- 
termined, p(r min )/pcrit = 9.4 x 10 5 , and p(r min ) / p crit ~ 
1.6 x 10 6 . These values can be combined with the require- 
ment of mass conservation to place an upper limit to the 
inner asymptotic slope of the density profile, f3 < 3(1 — 
p(r m in)/p(?"min)) = 1-2, although it is possible that the slope 
may actually become even shallower near the centre, as sug- 
gested recently by Taylor & Navarro (2001). 

Our results thus argue against the very steep values for 
the asymptotic central slope (/? w 1.5) claimed recently by 
Moore et al. (1998, 1999), Ghigna et al. (1998, 2000), and 
Fukushige and Makino (1997, 2001). The reasons for this dis- 
agreement are unclear at this point, since there are substan- 
tial differences in the halo mass, numerical techniques, and 
cosmological model adopted, which hinder a direct compar- 
ison between our results and theirs. For example, the work 
of Moore et al. (1998) and Ghigna et al. (2000) differs from 
ours in mass scale (they simulated a galaxy cluster while we 
target a galaxy-sized halo) and in cosmology (they adopted 
an Einstein-de Sitter CDM cosmogony, whereas we adopt 
the ACDM model). 

Finally, the difference between the conclusions from var- 
ious authors may just reflect the fact that each group ap- 
plies different criteria to the identification of the regions 
deemed trustworthy. We note that models with the very 
steep (/3 ~ 1.5) inner slopes proposed by the Moore et al 
group and with the shallower slopes that we find here are 
almost indistinguishable if we restrict our analysis to radii 
> 2% of the virial radius. Probing radii within the inner 1% 
of the virial radius seems required to shed light on this con- 
troversy. Further simulation work with resolution adequate 
to address this issue in detail is currently underway. 
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APPENDIX A: THE GENERATION OF 
COSMOLOGICAL INITIAL CONDITIONS 

Periodic boundary conditions are usually adopted in cosmo- 
logical simulations for reasons of convenience. The assump- 
tion of periodicity implies that the simulation volume as a 
whole has to have precisely the mean density, a requirement 
that places restrictions on the size of the region and on the 
redshifts at which a particular simulation may be consid- 
ered reliable. On the other hand, with periodic boundaries 
the density field can be expanded as a sum over a discrete 
set of periodic plane waves. For a simulation volume which 
is cubic, the Fourier transform of the density field has the 
form of a cubic grid in Fourier space. The discrete nature of 
the power spectrum thus makes it easy to set up Gaussian 
density fields. 

The aim of our initial conditions generating procedure 
is to provide a particle realization of a Gaussian density field 
with the chosen power spectrum, P(k), on scales and at red- 
shifts where linear theory is applicable. Our procedure fol- 
lows closely that described in Efstathiou et al. (1985), where 
further details may be found. As in Efstathiou et al. (1985), 
we use the Zel'dovich approximation to perturb particles 
from a uniform cubic grid, 



x(t) 



&(*W(q) 



(Al) 



where x is the comoving Eulerian coordinate of the parti- 
cle, q is the Lagrangian coordinate denoting the particle's 
unperturbed position in the grid, b(t) is the linear growth 



factor, and rp is a function that describes the spatial struc- 
ture of the density field. The function ip can be expressed in 
terms of the acceleration field at time t, 



<Mq) = - 



F(q,t) 



ma 2 (ab + 26a) ' 



(A2) 



where F is the force field, a(t) is the expansion factor, m is 
the particle mass, and a dot denotes a time derivative. 

In practice, a realization of the desired fluctuation dis- 
tribution is created in Fourier space, with random phases 
and normally distributed amplitudes for the real and imag- 
inary components of each mode. We then multiply by an 
appropriate Green's function, and transform back to obtain 
the potential on a spatial mesh. This potential is differenced 
to obtain F(q, t) which, together with eq. Al and eq. A2, 
gives the displacement field required to generate the desired 
density fluctuations from a uniform distribution. 

Once the displacements from the unperturbed positions 
have been computed, velocities are assigned to the particles 
assuming that only growing modes are present. The pecu- 
liar velocity is then simply proportional to the displacement 
vector, 



x = -&v(q) 



(A3) 



In cases such as CDM, where there is significant power on 
all scales, it is important to avoid unrealistically large initial 
velocities that may result from large amplitude fluctuations 
on small scales. Thus we assign peculiar velocities only after 
recalculating the accelerations using the perturbed particle 
positions and using eq. A2 to re-estimate V>(q). 

We use a cubic grid distribution of particles to represent 
a uniform density distribution for all simulations reported 
here, but it is also possible to use a 'glass' for the unper- 
turbed configuration. As discussed by White (1994), this 
is a better choice for highly aspherical perturbations, and 
avoids artifacts that arise from the existence of 'preferred' 
(Cartesian) directions in cubic grids. This is especially im- 
portant when attempting to simulate very low mass halos 
in CDM cosmogonies, since on those scales the mass fluc- 
tuation spectrum is nearly 'flat' (P{k) approaches fc ) and 
collapse proceeds almost simultaneously on many different 
mass scales in a network of sheets and filamentary struc- 
tures. 



APPENDIX B: 
TECHNIQUE 



MASS REFINEMENT 



As discussed above, simulations of periodic boxes are only 
reliable provided that the box is large enough so that per- 
turbations on scales comparable to the box size are still 
linear by the present time. This sets a minimum size for 
periodic boxes designed to be run to z = in a ACDM 
universe. For example, in the case considered in this paper 
(see § 2.2) the size of the periodic box is Lb ox = 32.5 ft -1 
Mpc (Mbox = 9.533 x 10 15 Q h' 1 M ), and the variance at 
z = is already <T(M bo x) ~ 0.3, at the limit of what may be 
used to obtain a good representation of large scale structure 
in this cosmogony. Clearly boxes smaller than 32.5 h^ 1 Mpc 
cannot capture the correct statistical properties of the dark 
matter distribution at z = 0. Our original low-resolution 
simulation was carried out with 128 3 particles. 



© 2002 RAS, MNRAS 000, 1-19 



Inner halo structure I: convergence 23 



Even if 512 3 particles were used in such a box (at 
the limit of what is possible with today's largest super- 
computer if many thousands of timesteps are needed) the 
mass per particle in a 32.5 ft -1 Mpc box would be m p = 
7.1 x 10 7 f2o/i _1 Mo, and a galaxy-sized, 10 12 h~ 1 M Q halo 
would only contain slightly more than 10, 000 particles. A 
dwarf galaxy halo would have fewer than 1,000 particles. 
Clearly, a different technique is required in order to improve 
the mass and spatial resolution of the calculation while at 
the same time accounting properly for the effects of large 
scale structure. 

The technique most widely adopted so far selects a few 
systems identified from the final configuration of the periodic 
box and resimulates the whole box, with coarser resolution 
everywhere except in the selected regions. This technique 
has been used in a number of cosmological simulations (see, 
e.g., Katz & White 1993, Navarro & White 1994, Evrard, 
Summers & Davis 1994, Moore et al. 1998), and has be- 
come common in high-resolution simulation work targeted 
at individual systems. The price one pays with this pro- 
cedure is that to build a statistically significant sample of 
halos entails running many different simulations and there is 
always the possibility of introducing biases during the selec- 
tion procedure. Having identified a halo in the periodic box 
for resimulation, all particles within ~ 2 r2oo from its centre 
are traced back to the initial conditions and their positions 
on the original cubic grid are recorded. A box of size Lsbox 
enclosing all of these particles is then defined. 

A displacement field is generated for iV s b x = 256 3 par- 
ticles in this new box using a two-step procedure that allows 
for inclusion of fluctuations from the original periodic box. 
In the first step, displacements for the iV s box particles are 
calculated using the same Fourier representation as in the 
original box, except for the contribution from wavelengths 
shorter than a characteristic scale, di cu t- Typically, di cut is 
chosen to be the shortest wavelength in the original box, 
dicut = 2L box /A b 1 ^ = 2 x 32.5/128 /i _1 Mpc ~ 0.5 ft" 1 Mpc, 
which is the Nyquist wavelength of the low-resolution par- 
ticle grid. We truncate the waves at a boundary which is 
cubical in Fourier space. 

It is important to ensure that the displacements due to 
the long wavelength Fourier components are applied to the 
high resolution particles in a sufficiently smooth fashion to 
avoid introducing significant spurious power. Computing the 
displacements by simple finite differencing of the potential, 
as is the case for the large periodic simulation box, is inad- 
equate in this context unless an impractically large mesh is 
deployed. A better way is to compute the individual compo- 
nents of the displacement field one at time, using the appro- 
priate Green's functions, and to interpolate the displacement 
components themselves, by trilinear interpolation to the in- 
dividual particle positions. The use of trilinear interpolation 
ensures that the displacement field is continuous-which in 
itself is sufficient to avoid spurious non-linear features being 
introduced. The larger the mesh used the more accurate is 
the interpolation. For the simulations reported here a 512 3 
mesh was used and proved satisfactory. 

In the second step, fluctuations are generated on scales 
smaller than di cu t, down to the Nyquist frequency of the 
high-resolution box. The new displacement field is periodic 
within Z/ s box, and can be vector added to the large-box dis- 
placements in order to obtain final perturbed positions for 



all particles within the high-resolution box. Trilinear inter- 
polation is once again used to assign the short wave compo- 
nents of the displacement field to the particles. Peculiar ve- 
locities proportional to the displacements are then assigned 
using the Zel'dovich approximation and assuming that only 
growing modes are present. 

Following this procedure, a realization of the displace- 
ment field of 256 3 particles is created and stored for each 
halo. Finally, the high-resolution box is inserted in the large 
periodic box after removal of all overlapping particles Not 
all particles in the small box will end up near the system of 
interest, so the location on the original grid of selected parti- 
cles is used to identify an 'amoeba-shaped' region within the 
cube that is retained at full resolution. Regions exterior to 
the 'amoeba' are coarse-sampled using particle masses which 
increase with distance from the region of interest (Figure 2). 
The sampling is typically done by binning together cubes of 
2 3n neighboring particles from the initial grid (where n is an 
integer). This allows us to concentrate numerical resources 
within our selected object without compromising the con- 
tribution from larger scales to the tidal field acting on the 
system. Because of the non-spherical nature of the collapse 
of dark halos, accurate simulation of the formation of a sin- 
gle system incurs a significant overhead. Even after all this 
optimization, at most 1 in 3 particles in the amoeba region 
ends up within the virial radius of the system considered. 

The success of our procedure may be gauged by com- 
puting the power spectrum from the displaced particle posi- 
tions and comparing it with the theoretical power spectrum 
that we are trying to generate. Figure 1 shows the desired 
theoretical power spectrum, the power spectrum measured 
from the parent simulation, and the power spectrum mea- 
sured from a high-resolution box created in the manner out- 
lined above. The power spectra are shown at z = 49. In 
this case, the high resolution box is 5.08 Mpc on a side 
and is sampled with 256 3 particles, with individual masses 
of 6.5 x 10 5 h~ Mq. The power spectrum of the small box 
is actually determined for a cube of 4.3 h Mpc excised 
from the middle of the high resolution region. The excised 
region would contain 216 3 particles if it had precisely mean 
density. The density field is assigned to a 432 3 mesh using 
a cloud-in-ccll assignment scheme and periodic boundary 
conditions forced. Forcing periodicity does not significantly 
distort the power spectrum for modes small compared to 
the fundamental mode of the cube. The power spectrum is 
then computed from the Fourier transformed density field. 
The power from individual modes is binned in shells of con- 
stant cubical wavenumber (k cubica i = max.(\k x \,\k v \,\k z \)). 
Plotting the power spectrum using the cubical wavenumber 
highlights discrepancies more sharply than the more usual 
spherical binning. The good agreement between the theoret- 
ical power spectrum and that measured in our realizations 
gives us confidence that our simulations faithfully follow the 
formation of a dark matter halo in the ACDM cosmogony. 
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